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ABSTRACT 


This work develops a simple, second“order paramuterization of the water 
fluxes at a landsurface for use as the appropriate boundary condition in gen- 
eral circulation models of the global atmosphere. The derived parameteriza- 
tion incorporates the high non-linearities in the relationship between the 
near-surface soil moisture and the evaporation, runoff and percolation fluxes. 

Based on the one-dimensional statistical-dynamic derivation of the annual 
water balance developed by Eagleson (1978), it makes the transition to short- 
term prediction of the moisture fluxes, through a Taylor expansion around the 
average annual soil moisture. 

A comparison of the suggested parameterization is made with other exist- 
ing techniques and available measurements. 

A thermodynamic coupling is applied in order to obtain estimations of the 
surface ground temperature. 
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NOTATION 


Symbola 

A 

A 


H 


w 


ij 


d 

d. 


Definition and Dlmenolons 

albedo [-] 

gravitational infiltration rate as modified by capillary rise 


from water table, 
biomass production 

coefficient for the sensible heat transfer, 
coefficient for the water vapor transfer 
specific heat of water vapor at constant pressure 
pore disconnectedness index 
specific heat of soll-w<? u '.r system 

moisture transfer coefficient between layers i and j 
depth to which the diurnal moisture cycle extends 


[LT"^] 

-2 -1 
[ML T ] 

[-] 

r-] 

[L^T'^deg"^] 

[-] 

[hV^deg"^] 

[LT"^] 

[L] 


a soil depth influenced by the diurnal soil-moisture cycle [L] 


dlffuslvity index of soil 

a soil depth influenced by the annual temperature cycle 

drying power of the air 

annual potential evapo transpiration 

annual actual evapotranspiratlon 

annual storm surface retention 
exfiltration parameter 

average annual actual evapotranspiratlon rate 
average annual potential evapotranspiratlon rate 
actual evapotranspiratlon rate 
potential evapotranspiratlon rate 
transpiration rate from vegetation 


[-] 

[L] 

[FL“V^] 

[L] 

[L] 

[L] 

[-] 

[LT~^] 

[lt“^] 

[LT"^] 

- 1 , 


[LT 
[LT“^] 


10 


e soil evaporatxon late 

s 

saturation vapor pressure at surface temperature 

e vapor pressure of the air at screen height 

o 

exfiltration capacity of soil 

f^ Infiltration capacity of soil 

G heat flux into the soil 

G gravitational infiltration parameter 

H sensible heat 

H sensible heat 

s 

I infiltration on soil surface 

c 

infiltration under vegetated surfaces 

i rainfall rate 

J evapotranspiration efficiency 

K(l) saturated hydraulic conductivity 

hydraulic conductivity between layers i and j 

atmospheric heat conductance 

K surface heat conductance 

s 

k(l) saturated intrinsic permeability 

plant coefficient 
k^ soil thermal dlffuslvity 

L latent heat of vaporization 

L Monln-Obukhov length 

M vegetal canopy density 

equilibrium vegetal canopy density 
m^ rainy season length 


[FL"^] 

[LT“^] 

[FL'^] 


[-] 

[FL~^] 

[FL"^] 

[-] 

[LT"^] 

[FT"^deg“^] 

[FT“^deg“^] 


[-] 


[L^T-^ 


1 

] 


[L] 

[-] 


t-] 

[T] 
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mean time between storms 

[T] 

nip 

mean annual precipitation 

[L] 




“v 

mean number of storms per year 

[-] 


mean storm intensity 

[LT”^] 

“h 

mean storm depth 

[L] 

m 

pore size distribution index of soil 

[-] 

n 

e 

effective porosity of the soil 

I-] 

n 

effective porosity of the soil 

[-] 

"a 

annual precipitation 

[L] 

p 

mean storm intensity 

« — 1 
1 

H 

p 

a 

a*'mospheric pressure 

CM 

1 

q* 

saturated atmospheric specific humidity 

[-] 


specific humidity of the atmosphere at screen 

elevation [-] 


soil moisture flux between layer i and j 

[LT"^] 

^Vb 

bulk Richardson number 

[-] 

R 

®A 

annual groundwater runoff 

[L] 

h 

annual surface runoff 

[L] 

^A 



R 

n 

tiet radiation at the surface 


R 

gas constant 

2 -2 

[L T deg 

r 

a 

atmospheric diffusion resistance 



surface diffusion resistance 


S 

e 

exfiltration "desorptivity" 



infiltration "sorptivity" 


s 

s 

relative humidity at the evaporating surface 

[-] 

S 

g 

saturation ratio of the surface atmosphere 

[-] 
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u 

u 

a 

w 

W 

g 

w 

n 

Y 


max 


y 

Z 


average soil moisture at the surface layer 
average annual soil moisture at the surface layer 
soil moisture concentration at time k 
critical value of soil moisture 


[-3 

[-] 

[-3 

[-3 


average annual atmospheric temperature [deg3 

air temperature at screen height [deg] 

one year [T] 

ground temperature at the surface [deg] 

surface temperature [deg] 

soil temperature at the depth of the vapor source [deg] 

mean soil temperature of layer of depth d 2 [deg] 

time when the surface becomes saturated after a precipitation [T] 

time when the surface becomes dry during an evaporation period [T] 
storm duration [T] 

time between storms [T] 

moisture uptake by plants 
wind speed 

upward capillary rise velocity from the watc-j table 
surface soil moisture [~] 

maximum value of soil moisture for which surface runoff is equal to zero [• 
total yield [L] 

percolation rate 


[L^^] 




average annual percolation rate 
surface runoff rate 

average annual surface runoff rate 
toTi.^l yield rate 
surface layer thickness 




[L] 
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z 


a 

6 

Y 

A 

6 


n 

0 



K 


X 


a 


w 


<*>1 

g 

H' 

n 


distance above the soil surface 
screen height 
surface roughness 

reciprocal of average storm Intensity 
reciprocal of average time between storms t^ 
the psychrometrlc constant 

slope of the saturation vapor pressure-temperature curve 

reciprocal of average storm duration t^ 

reciprocal ot mean storm depth m^^ 

volumetric moisture content 

field capacity 


[L] 

[L] 

[L] 


[rt] 

[t’^] 


[FL“^deg"^] 

[FL“^deg“^] 

[t“^] 


[-] 


[-] 


shape factor of Gamraa-dlstrlbuted rainstorm depths 

parameter of Gamma-dlstrlbuted storm depth 

potential humidity 

mass density of air 

mass density of water 

density of soil-water system 

capillary Infiltration parameter 

one day 

dimensionless desorption dlffuslvlty of soil 

dimensionless sorption dlffuslvlty of soil 

soil matrix head 

soil matrix head 

groundwater recharge potential 


[-] 

M 

[-] 

tFL~S^] 
[FL'S^] 
[ ] 

[T] 

[-] 

[-] 

[L] 

[L] 

[-] 
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CHAPTER 1 


Introduction 

1.1 BactcRround 

Current global atmospheric general circulation models use very complex 
numerical techniques to solve the hydrodynamic and thermodynamic equations 
of motion in the atmosphere, but they generally treat the land surface ther- 
mal and moisture boundary conditions in a rather simplistic way. This study 
attempts to provide an Improved land surface boundary condition that increases 
the physical fidelity while maintaining computational practicality. 

There are many difficulties involved in such a parameterization. First, 
there are problems related to inhomogeneities of the system's inputs, such as 
precipitation and of the system parameters, such as soil properties. Because 
of the great difficulty in defining the interactions of the microscale and 
the macroscale dynamics and representing these in a computationally efficient 
way, the problem of spatial variability is treated by considering a lumped one- 
dinensional system, having effective areal parameters. A second problem is 
that of formulating the appropriate differential equations, which will account 
for the exchange of water and heat between the soil and the atmosphere. Spe- 
cial problems that arise here are those concerning the time scales of the phys- 
ical processes, the number of parameters used and the selection of conceptual 
models representative of the real processes. 

1.2 Objectives 

The objectives of this work are; 

1. To derive analytical expressions for the evaporation and yield rates 

which can be applied in dynamic mass balance equations for short term 

prediction of the soil moisture in the root zone. 
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2. To minimize the number of parameters necessary to implement this 
parameterization and to determine the inputs and observations required 
to operate the model. 

3. To compare the results with those obtained from other models, which 
use either detailed numerical techniques or different types of paramet- 
erization. 

4. To perform sensitivity analyses with respect to the critical soil 
parameters. 

5. To estimate the ground surface temperature. 

It must be noted, that the presence of snow is not taken into account in 
this research. 
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CHAPTER 2 


Literature Review 
2.1 Landsurface Parameterization 

According to Eagleson [1981], there are six basic elements of the hydro- 
thermal system at the earth surface, which must be parameterized; 

1. The rate of potential evapotranspiration e , which is a function of 

P 

the incoming short and long wave radiation to the system, the wind speed, 
the surface roughness, the vapor pressure of the air, the temperature of 
the air and of the ground. The dependence of e^ on the ground tempera- 
ture generates a feedback between the soil and the atmosphere, and be- 
comes a major coupling factor between the heat and moisture fluxes, Be- 
cauf e of the creation of an atmospheric boundary layer due to the air 

flow close to the surface, e becomes also a function of the extent of 

P 

the upwind evaporating surface. 

2. The actual evapotranspiration rate e^, which is a function of the 
available soil moisture, the soil properties and the vegetation cover. 
The value of e^ is limited from above by the value of e^, i.e., the cap- 
acity of the atmosphere to remove vapor from the surface. The following 

general expression relates e with e : 

P 

j = ^ = f(s^ e , t; soil and vegetation), J ±l 

% p 

where J is called the ’’evaporation efficiency”, 

3. The water yield rate y, which is divided into two components, the 
surface runoff rate y and the percolation to the water table y . The 

S o 

yield rate is functionally related to the following parameters: 
y = f(s, precipitation input, t; soil properties and storage capacity 
of the surface layer) . 


17 



4. The surface temperature T , which is dynamically related to the net 

S 

radiation R , the latent and sensible heat losses from the soil, and 
n* 

the heat storage capacity of the surface soil layer. 

5. The surface moisture retention capacity e^, which can become im- 
portant for certain types and density of vegetation cover and also under 
very arid conditions. 

6. The soil moisture layer thickness Z^, which consists of the portion 
of the soil close to the surface, where changes in moisture and heat con- 
tent are concentrated. This zone is usually taken equal to one meter, 
but in fact it should be determined by the root-zone depth and by the 
soil and climatic properties of the area under investigation. 

An overview of the methodologies proposed by prior investigators to mod- 
el the above elements, will now be made following that of Eagleson (1981). 


1. Potential Evapotransplration rate ep 

The concept of potential evapotransplration, first introduced by Thornthwaite 
[1948], refers to the capacity of the atmosphere to remove vapor, under given 
meteorological conditions, when there is unlimited water supply from the soil 
and the ground surface is wet. 

The basis of the recent approaches for estimating e^ is Penman’s [1948] 
equation in the modified form: 


where 


E 

a 


and 


A(R -G) E (x) 
D T e = —2— + v-S 

^w *p A + y 'A + Y 


Q L 


r 


a 


[q*(T ) - q ] 
a a 


( 2 . 1 ) 
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R = net radiation near the surface, 

n 

G = net heat flux into the ground. 

A = (de*/dT) , the slope of the saturation vapor pressure-temperature 

curve . 

Y = the psychrometric constant. 

L = latent heat of vaporization 

E = the "drying power" of the air 

q*(T^) = saturated atmospheric specific humidity at air temperature 

q = specific humidity of the atmosphere at screen elevation 

r^ = atmospheric diffusion resistance. 

Equation (2.1) was applied in many theoretical and experimental studies 
by investigators such as Penman and Schofield, 1951; Penman, 1956; Tanner and 
Petton, 1960; Slatyer and Mcllroy, 1961; Monteith, 1965, 1973; Rljtema, 1965; 
Van Bavel, 1966; Kohler and Parmele, 1967; Thom and Oliver, 1977. 

The two- term structure of Equation (2.1) helps to point out the influence 
of large-scale advection. The first term corresponds to a lower limit to evap- 
oration from a moist surface under steady-state conditions. When such con- 
ditions have been established, the value of q tends to reach the saturation 

a 

value. This can happen in the downwind direction of a wet surface of infinite 
extent, where evaporated moisture will be advected. The second term of Equ- 
ation (2.1) represents the drying power of the air E (x) . It takes it’s 
maximum value at the beginning of the uniform surface (x=0) , where the air is 
dryest. It was found by McNaughton [1976], that E^(x) decreases exponen- 
tially with distance x. 
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2. The Actual Evaporation Rate st 


When the limiting factor for evaporation is not the available energy 
supplied to the system, but is the amount of soil moisture within the sur- 
face layer, then evaporation control passes to the soil, in this case, the 
evapotranspiration rate becomes dependent on the value of soil moisture, the 
soil properties and the vegetation cover, which together influence the capa- 
city of the soil to deliver moisture upwards. Thus we can write; 

e = f(e , a; soil, vegetation). The methods developed in order to deter- 
P 

mine this function can be categorized as follows; 

a. Empirical parameterizations 

Those involve long-term average relationships between precipitation 
and evaporation, which are derived from simple water balance equations 
applied to various catchments, by equating the total streamflow at the 
end of the catchment to the total yield produced in the basin. Such em- 
pirical relations are of no importance if one is interested in under- 
standing the dynamics that govern the physical process, the Interaction 
between the system's elements and their response to different hydrologi- 
cal and atmospheric conditions. As references one can mention the works 
by Schreiber [1904], Ol'dekop [1911], Tara [1954], Budyko [1956], Pike 
[1964], Budyko [1971]. 

b. Divergence of Atmospheric Vap^r Flux 

Rasmussen [1977]; derived a steady-state long-term regional atmos- 
pheric water balance over an area in space, by considering horizontal 
water vapor fluxes measured with probes well distributed in space. He 
used surface precipitation observations to close the water balance and 

estimate the long-term spatially-averaged actual evapotranspiration. 
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c. Advection-Aridity Approach 

Brutsaert and Strieker [1979] assumed that the second term of Equa- 
tion (2,1) represents the effect of larger-scale advection. They also 
used the concept of symmetry between potential and actual evapotranspir- 
ation introduced by Bouchet [1963] and the corresponding value of e^ 
for conditions of minimal advection suggested by Priestley and Taylor 
[1972], to derive: 


a(R^ - G) 


R - G + +E 
n A a 


-1 


( 2 . 2 ) 


where a ^ 1.26 

This methodology is called the "advectlon-aridity" approach. The 

time scale in Equation (2.2) is arbitrary, although they found it to 

work satisfactorily for dally values. It must be noted that difficulties 

are encountered in estimating the wind function which enters into the 

calculation of E and also advection effects were assumed to be gener- 
a 

ated only due to the regional shortage of moisture supply at the surface, 
d. Soil Moisturj Surrogate 

Attempts have been made to derive equations for the actual evapo- 
transpiration rate, by introducing a soil-moisture related surface par- 
ameter. All equations of this type are of the general form: 


i + .y B 

A+y 


(2.3) 


where 


i. B = K /K , Slatyer and Mcllroy [1961] 

Sl S 

K^,K^ = atmospheric and surface heat conductances. 
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ii. B >= Monteith [1965] 

r = atmosiiheric diffusion resistance 
a 

r “ surface diffusion resistance, which was related to the 
c 

evaporation rate and to the difference between the vapor pres- 
sure at the leaf surface and it's saturated value at leaf sur- 
face temperature. 

iii. B = (1-S )/S , Barton [1979] 

S 8 

where S = relative humidity at the evaporatin'* surface, which 
s 

for non saturated surfaces can be very different from the rel- 
ative humidity at the ground surface. 

A special difficulty encountered in expressions such as Equation 
(2.3) is that they are not designed to represent effective areally aver- 
aged values for the atmospheric temperature and the value of b. 

Tanner and Fuchs [1968] derived a more general equation for e^, 
which does not assume any particular diffusion model for the leaf or 
other surface and does not include internal resistance to heat diffusion. 
The model they used is: 


Ep- [ pCpA/ (A+ y) ] (T^-T^) / r^ 
^ " l+[Y/(A+Y)](r^/r^) 

where 


(2.4) 


r = internal resistance of the soil surface layer to the transport 
of water vapor [sec.m ] 

T^ = surface temperature [*K] 

T^ = soil temperature at the boundary between the dry and moist layer 
in the soil, where the vapor source is. 


r 
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There are also studies which derived expressions for the transpiration 
from plants. A category of these assumes knowledge of plant physiology, which 
can reduce their applicability for macroscale parameterizations, due to their 
inability to capture dynamic interelationships among the various spatially 
variable system components. Studies of this type include those by Van de Honert 
[3948], Cowan [1965], Rijtema [1965], and Federer [1977], 

Assuming the same albedos from the vegetation and the wet bare soil, 
Shuttleworth [1979] proposed a relationship similar to Equation (2.3), where 
he replaced e^ by transpiration rate from the plant. 

Eagleson [1981] applied ecological optimality hypotheses to water-limited 
Natural Soil-Vegetation systems, to derive the following equation for the 
average evapotranspiration efficiency; 



' 0.11 + 2.22M - 1.87M^ + 0.54M^ 

o o o 

0.11 + 1.25M + 0.27M^ - 0.63M^ 

0 0 0 


, k “ 1 

’ V 


9 


k 

V 


n 


0.7 


(2.5) 


where 


= percentage of vegetation cover 

. ^Tv 

and k “ -=r— = plant coefficient. 

V 6p 

e. Moisture Accounting Models 

Most of the atmospheric general circulation models currently in use, ap- 
ply a surface boundary conditions which incorporates an evaporation-soil- 
moisture relationship of the linear Thornthwaite-Budyko type. This relation 
has the form: 



0 > 

0 1 0 1 0fc 


( 2 . 6 ) 
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where t)f^ is the soil field capacity, i.e., the upper limit of soil 
moisture for which water can be stored in the soil without drainage due to 
gravity. 

A major difficulty encountered in formulations such as Equation (2.6) is 
the definition of the field capacity since there is no rigorous justi- 

fication for the correct value of this parameter. 

As it is pointed out by other investigators (Philip, 1957; Hlllel, 1971; 
Lowry, 1959), the relation between evapotranspiration and soil moisture can 
be highly non-linear, due to the influence of soil properties and vegetation 
cover, which play a dominant role during the exfiltration process. This non- 
linear relation was also theoretically supported and verified by Eagleson 
[1978], and it is the one used for the purposes of the present study. 

Eagleson [1978] used the probability distributions of the Independent 
climatic variables to obtain the derived distributions of the dependent ele- 
ments of the water-balance, through a physically-based model of the natural 

process. With this statistical-dynamic approach he derived an expression of 

Bt 

the long-term annual average evapotranspiration efficiency ■=— as a function 
of the long-term averages of soil moisture and climatic parameters, soil 
properties and vegetation cover. This relationship will be presented with 
more details in Chapter 3. 

3. Water Yield Rate y 

There are empirical equations which relate the long-term average annual 

yield with annual precipitation and potential evapoi'atlon (Lettau, 1969). 

Another way of estimating the total long-term yield is by equating it 

with the long-term time integral of the streamflow in the catchment. This 

assumption can lead to errors especially under very arid conditions, where 

the groundwater flow does not appear as streamflow. 
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Empirical raodels developed by Budyko [1971], and Arakawa [1972], osti“ 
mate the short-term yield as a function of soil moisture, precipitation, po- 
tential evapotranspiratlon, soil porosity, and surface layer thickness Z^. 


4. Surface Temperature 

The ground temperature T is an Important parameter for determining 
sensible and latent heat fluxes from the soil to the atmosphere. Complicated 
numerical models exist, which solve the coupled moisture and heat transport 
equations in porous media. Recent studies include the work by Philip and 
DeVries [1957], Sasamori [1971], Rosema [1975], Benoit [1976], and Milly 
[1980]. 


For application in GCM's, more simplified methods for estimating T are 

S 

needed. The most commonly used among those methodrlogies as presented by 
Deardorff [1978] are: 

a. Insulated surface (Oates et al. [1971]; Manabe et a ''74]. The 
heat flux into the surface G is taken equal to zero and .. ' bal- 

ance equation at the surface l.e., 


R (T ,t) - H(T ,e) - p L e„(T ,t) G 
n g’ 8 w T g* 


(2.7) 


must be solved for T , given that the other elements of the equation 
are known. 

b. Dependence of G on the sensible heat H. Kasahara and Washington 
[1971] assumed G = 4 h and solved Equation (2.7) for T . 

j 8 

c. Dependence of G on the net radiation R^. Nickerson and Smiley 

[1975] assumed G = -0.19R , when R <0 (down) and G “ -0.32R , when 

n’ n n’ 

R^>0 (up) and solved Equation (2.7) for T^. 

d. Bottom- insulated single soil layer. Arakawa [1972], Corby [1972], 


and Rowntree [1975] applied the following equation; 
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whore 


^8 " 


d 


1 


a 


and dj^ ° 
where 


p^(0) ■=> density of soil-water system 
c (0) => specific heat of soil-water system 

depth of the soil layer influenced by the diurnal temperature 
cycle 


° soil thermal diffusivity 
° one day 

e. Fr"*ce-re8tore method. Bhumralkar [1975] and Blackadar [1976] 
developed a formula that contains the deep soil temperature T 2 of the 
following form; 

a G/(p c d.) - 2 tt(T - T,)/t. (2.9) 

dt S S 1 g 1 

where 

a t= I + — and T = T(6,t), 6=lcm 
d 6 

1 

Deardorff [1978] considered the case where (S—>-0 and Lin [1980] 
considered the 6 layer thickness effect to be weaker than that of 
Bhumralkar [1975] by setting; a=l + r 

Tests performed by Deardorff [1978] and Lin [1980], proved that the 
force-restore method gives reasonably accurate results for estimating 
the ground surface temperature. 

Sasamori [1970], developed a model by using the thermodynamic equil- 


ibrium relation; 
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( 2 . 10 ) 


° exp[gi|) (0)/RT ] 

o b S 

wher*^ 

S a saturation ratio of the surface atmosphere 

(0) 13 soil matrix head as a function of 0 

h 

R “ gas constant 
That way he provided a coupling between the energy and mass conservation 
equations and the local thermodynamic equilibrium of temperature and 
humidity. Equation (2.10) could be solved for T in the case of soil con- 
trolled evaporation given thac the other terns are kno\m. A special dif- 
) Iculty could occur when the surface is saturated. 

Then T should be approximated by the temperature above the evapora- 
ting surface. 

5. Surface retention capacit y 

There Is a volume of precipitation moisture which is retained at the sur- 
face due primarily to the surface texture. There are empirical relations 
for estimating that capacity and a collection of them is given by Wigham 
[1975], Blake [1975], and others. It should be noted that this water loss 
must depend also on the amount of precipitation, it^s intensity and the dura- 
tion of interstorra periods. 

6. Thickness of Soil Moisture Layer Z,, 

This layer represents the depth from the surface within which big changes 
in soil moisture and heat content can occur due to forcing from the atmosphere. 
By using Philip’s [1969] infiltration theory, Eagleson [1978] showed that 
this capillary penetration depth is of the order of one meters Clearly, this 
depth would be a function of the soil-type, the timing of precipitation events 
and the depth of the root-zone system. 


The diurnal thermal penetration depth has been found to be of the same 
order of magnitude. 

Budyko assigned a value of ° Im and Arakawa assumed nZ^ ° 10cm. Gates 
et al. [1977] suggested = 30 cm and Shukla suggested 0fj,Zj- ® 10cm. 

Clearly, a rigorous justification of the appropriate value of Zj. does 
not exist and it's value is chosen rather arbitrarily. A sensitivity analysis 
is needed, in order to define the critical parameters that Influence Z^. 

2.2 Water Balance Models 

The existing water balance models can be divided into three categories: 

a. Empirical (Thornthwaite and Mather, [1955]; Lettau, [1969]). 

b. Phenomenological ("The Stanford Watershed Model", Crawford and 

Linsley [1966]; Holtan et al. [197A]; Peck [1976]). 

c. Dynamic (Eagleson [1978]). 

For short-term predictions of soil moisture the Deardorff [1978] and 
Lin [1979] models are mentioned here, since results from their methodologies 
are compared with those obtained in this work. 

Deardorff [1978] applied the "Force-restore" method to determine surface 
moisture and temperature. He used a value of the bulk soil moisture In the 
upper half -meter of the soil and wrote an equation of the form: 



(e -i) (W -W.) 

-n ^ _ r K D 

^'1 p d/ '"2 T 
w 1 


0 < W, 
— b 


< W 


max 


( 2 . 11 ) 


where 


28 


1 


w 

max 


1 = storm intensity 

W = surface soil moisture 

g 

C, = f(W /W ) 

1 ' g max 

= constant 
T “1 day 

depth to which the diurnal moisture cycle extends. 

maximum value of soil moisture for which surface runoff is equal 
to zero. 

Transpiration from vegetation was Included by using a generalization of the 
Monteith and Szeicz [1962] function for evaporation, which incorporates par- 
ameters related to plant physiology. The complicated representation of vege- 
tation can make this approach difficult to apply in an actual situation. 
Gravity is ignored and actual evaporation is derived from a Budyko type lin- 
ear relation. 

Lin [1979] developed a deteministic model to be dynamically coupled 
with a GCM and the groundwater zone. The ground is represented by a surface 
layer of 10cm depth, an intermediate 40cin layer and a deep layer which con- 
tains the groundwater zone. 

The equations describing the system’s dynamics are given by: 

Surface layer di : 


d^(dG^/dt) = (Ig-e^/p^) (1-M) + - q^2 

Intermediate layer d? : 

d2(d02/dt) = -U 2 M + q^^2 ~ ^123 

Deep layer : 

0^ = f(t), with very slow variations with time. 


( 2 . 12 ) 


(2.13) 


where 
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I ,I 
s’ c 


M 


Ui,U2 


‘ij 


= Infiltration on soil and under vegetated surfaces 
= soil evaporation rate 
=» percentage of vegetation cover 
“ moisture uptake by plants 
= moisture flux from layer i to j 


The potential evapotranspiration is derived from an aerodynamic equation 
and the actual evapotranspiration as a refinement of Budyko’s parameterization, 
by using the value of the field capacity, wilting point and some empirical 
constants. The surface retention capacity of the soil and vegetation is ig- 
nored. The infiltration capacity is estimated by applying Holtan’s [1974] 

B 

method, where s A0^ . The soil moisture flux > between adjacent layers 
is given by: 


where 


1, . = D, .(0.-G.) + K, . 
ij iJ J J ij 


(2.14) 


D 

K 


ij 

ij 


moisture transfer coefficients 
hydraulic conductivity of the soil 


From nu/aerical experiments, he derived reasonable results of hydrologic 
variables such as soil moisture, evapotranspiration and surface temperature, 
for various regions of the earth. He distinguished between the variables that 
can be obtained from remote sensing and which are 0^ , M, and T and those 
which cannot. The time and space varying parameters and create a 
special difficulty to implement in the model, since very often large scale 
measurements of those parameters do not exist. 
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CHAPTER 3 


Review of the Water Balance 

3.1 Introduction 

The theoretical background of this work is drawn from Eagleson's (1978 
a, b, c, d, e, f, g) water balance model. The physical model is one-dimensional 
and only vertical fluxes of water are considered in the soil column. The inputs 
to the system are of two types: 

a. Climatic variables, which are treated as independent random variables 

and are seven in total number. 

b. Soil properties, which are represented by three Independent parameters 

considered to be deterministic. 

The effect of vegetation is explicitly considered in the model through 
two parameters, the percentage of v jgetation cover M^, and the plant water use 
coefficient, k^. Vegetation is modeled to act as a uniformly distributed sink 
throughout the entire surface layer of thickness Z^, which continuously extracts 
moisture during the evapo transpiration period, at a rate regulated by the value 
of k . 

V 

The use of natural selection hypotheses and possible observations of the 
percentage of vegetation cover and total water yield from the basin, can signi- 
ficantly reduce the number of necessary input parameters to the model. Those 
techniques will be referenced with details in Chapter 5. 

Uncertainty in the model is introduced through the probability distribu- 
tions of several climatic variables. Precipitation events are simulated as in- 
dependent and identically distributed rectangular intensity pulses having 
Poisson distributed arrivals. The corresponding storm depth s h, are assumed 

as gamma distributed. The interstorm periods t, and storm durations t , are 
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taken to be exponentially distributed. The rate of potential evapo transpira- 
tion e^ and the air temperature T^, are set equal to their annual average 
values. The system dynamics for soil moisture movement are represented by 
Philip's (1969) infiltration and exfiltration equations. The averaged out- 
puts from the system, i.e., the actual evaporation E„, the surface runoff R 

i. s 

and the groundwater runoff R are calculated through the use of derived dis- 

s 

trlbutions. More details for the model's assumptions are described by Eagleson 
(1978 a). 

The mean annual water balance equation [Eagleson, 1978e] is given by: 


E(P^) {l - 
where 


E[ ] 



J 

E 

G 

a 

m 

T 

K(l) 

s 

o 

c 

T 

w 


e"G"2°r(a+i)a-°l = E[Ep *]J(E,M,K^) + m^K(l)s^'^ - Tw 
^ A 

» expectation operator 
= annual precipitation 

= annual potential evapotranspiration 

= annual actual evapotranspiration 

= /Ep * = evapotranspiration efficiency. 

= exfi'.tration parameter 
= gravitational infiltration parameter 
= capillary infiltration parameter 
= rainy season length 
= saturated hydraulic conductivity. 

® average annual soil moisture at the surface layer 
= pore disconnectedness index 
= 1 year 

= upward capillary rise velocity from the water table. 


(3.1) 


3,2 The Separate Elements of the Water Balance 


The basic elements of the water balance, which will be of use in the 
current landsurface parameterization are the following; 
a. Evapo transpiration 

Eagleson (1978d) derived the total evapotranspiration during an in- 
terstorra period by using an exfiltration analogy of Philip's [1969] in- 
filtration equation of the following form: 


where 


M = 


e 

V 

w 


S 

e 


where 


n 

e 


K(l) 


V(l) 


d 

m 


s 

o 


f - ^ S - M-e + W (3.2) 

e 2 e V 

vegetation fraction of surface 
vegetation transpiration rate 

velocity of capillary rise from the water table 

the exfiltration "desorptivity** which is defined as follows for a 
dry surface (s^=0) 

S = 2s 
e o 

= effective porosity of the soil 
= saturated effective hydraulic conductivity 
= saturated matrix potential of soil 
= dimensionless desorption diffusivity of soil 
= diffusivity index of soil 
= pore size distribution index of soil 

= initial soil moisture, constant throughout the surface boundary 
layer 


n^K(l)iJj(l)(t)^(d) 

Tim 


(3.3) 
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The boundary conditions associated with the exfiltration equation 
are interpreted as follows: 

In the beginning of the inters torm period, the surface soil moisture 
s will take a value s- == 1, so that the exfiltration rate f will be 
equal to the potential evaporation rate from a wet surface. We denote 
by the time it takes for the surface retention to evaporate and by t^ 
the time it takes for the surface to become completely dry. When 
t > t* + t^, then f^ < e^, anu f^ will be given by Equation (4.1). Evap- 
oration ceases at the time when = 0, or when a new storm begins. 
Transpiration from the vegetated surface assuming unstressed conditions, 
will take place at a constant rate e^, which will be givey by: = *^v^p’ 

where k^ « plant water use coefficient and e^ = bare soil potential evap- 
oration. 


t 

o 


The time t^, after which control passes to the soil was found to be 
equal to; 


2e ^(1+Mk -w/e ) 
P V p 


1 - M + 


I 


M^k + l(l-M)w/e 

__v E 

2(1+Mk -w/e ) 

V p j 


(3.4) 


Assuming exponential distribution of the time between storms t^ and 
a constant potential evapo transpiration rate e^ equal to it's annual av- 
erage value, Eagleson [1978b] derived the following expression for the 
average annual evapotranspiration efficiency J(E, M, K.^) s 
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OF poy:: 


II 

'p 


j = ^ = 1 

e, 


_ -j 

1_1 - M + Mk 1 

- [mU^ + (2C)^' 


1 + Mk^ + (2B)\ 


- (2E)^ 


Hl. H] 


(3.5) 


where 


= actual annual average rate of evapotranspiration 


B = 


1 - M 


M k + (l-M)w/e 
V _E 


- .2 


1 + Mky - w/e 2(1 + Mk^ - w/e^) 


C = I/(Mk^-w/ip)^ 

23n K(l)tKl)4>^ ^^2 

and E = ;; s 

— JL O 

TTine 

P 

where 3, Is the reciprocal of the average time between storms m . 

b 

b. Surface Runoff 

Assuming uniform intensity i(t) during a rainstorm Eagleson [1978e] 
applied Philip's [1969] infiltration equation, to represent the infiltra- 
tion rate by; 


^i “ I ^ ''o 


(3.5) 


for a saturated surface. 


and 


where 


A = ^ K(l)(l+sJ^) - w 

O L O 


S. = 2(l-s ){[5nK(l)i;j(l)$. (d,s )]/3mTi}'‘ 

1 O 10 


= infiltration sorptivity 


<t'^(d,s^)= dimensionless sorption diffusivity of soil. 


35 


When the surface is saturated (Sj^“l) , f becomes maximum and equal to the 

it 

infiltration capacity f^ . 

If the depth of surface retention is represented by h^, then infil- 
tration into the soil will start if t^ > h^/i. The initial infiltration 
rate will be equal to the stem intensity i. The continuous rise in the 
internal soil moisture will cause an increase of the surface moisture. 
When the surface becomes saturated, at time t + h /i, the inf ilt..ation 
rate will be given from Equation (3.6). Thus, for > t^ + h^/i, i > f 
and surface runoff is produced. It is found that: 


g2 ^ ^ 

*^o “ 2i(i-lT) 2'(i-A^~)j 

By solving the linearized diffusion equation with constant flux 

boundary condition, a similar expression can be obtained for t^ (Carslaw 

2 

and Jaeger, 1959), where the coefficient 1/2 is replaced by tt /16, 

By assuming the value of to be constant at it’s time- aver aged 
value, Eagleson [1978e] approximated the expected annual surface runoff 
R with the following function. 


E[R^ ] 

= e“^^'^r(0+l)/a'^ - E[E^]/m^ 
A 

where « „ 1/3 

■Snn K(l)ip(l)(l-s^) •I>^(d,s^)'‘ 

^ ~ 6lr6m 


G = [aK(l)/2] [1+s^] - aw 

E^ = storm surface retention 

ri = reciprocal of mean storm depth m^ 

6 = reciprocal of average storm duration t^ 

a = reciprocal of average storm intensity m 
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(3.8) 


Percolation to che Water Table 


C t 


The average annual groundwater runoff is given by [Eagleson, 1978f); 


where 


E[Rg ] = m^K(l) - Tw 


ra^ = mean rainy season length 
T = one year 


(3.9) 
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CHAPTER A 


The Short-Term Water Balance Model 

The purpose of this model is to make short-tem predictions of the soil 
moisture level within the surface layer, by taking into account the atmospheric 
boundary conditions at the surface. 

We assume that high moisture concentrntlons occur within a depth from 
the surface. That Implies that a portion of the infiltrated water during pre- 
cipitation will be stored within that layer of thickness Z^. Bel^w that depth 
Z^, water will percolate downwards due to gravitational forces. If the sur- 
face reaches saturation during the precipitation period, then runoff will be 
produced at a rate depending on che value of the internal moisture within the 
layer Z^. During evapotranspiration, the actual evaporation rate from the 
bare soil will be determined from the amount of available moisture within this 
layer and will be limited from above by the potential evapotranspiration rate 
Cp. Vegetation is assumed here to transpire at the potential rate under un- 
stressed conditions, which will be independent of the level of soil moisture 
within the layer Z^. 

We consider a vertical soil column in contact with the atmosphere and 

define as the moisture state variable s , the concentration of soil moisture 

0 

in the vegetation root zone Z^. In a one-dimensional representation, where 
only vertical fluxes are considered, the conservation of water mass equation 
can be written; 


where 


nZ 


ds 

dt 



( 4 . 1 ) 


e^ = evapotranspiration rate = f^^(s; climate, soil, vegetation) (4.2) 
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(4.3) 


y => yield rate ° f 2 (s; climate, soil) 

r a effective porosity of soil in vegetation root zone 

“ thickness of vegetation root zone (cm) 

i => rainfall rate (cm/sec) 

e^ “ evapotranspiration rate (cm/sec) 

y “ yield rate (surface runoff plus percolation below the root zone) 
(cm/ sec) 

Equations (4.2) and (4.3) establish a non-linear relation between e^, y, 
and s. In order to make the transition from the long-term time averaged 
values and y to their short-term values, e^ and y vary linearly around their 
long-term averages, with respect to the value of s. Thus, the non-linearities 
in the relations between e^, y, and s will be incorporated into the model 
through a Taylor expansion of those functions, around the annual average soil 
moisture s^. By performs that expansion, a transition is made from the long- 
term average values of those rates as they appear in the water balance derived 
by Eagleson [1978], to their short-term values. 

The water balance can be written in the following form, which is more con- 
venient for this purpose: 

1 = "J + e"^"^“r(0+l)o"‘^ = (4.4) 

o 

where 

H = potential humidity = E[P.]/E[E ] 

C “ groundwater recharge potential « m K(i)/E[P.] 

I A 

ep 

J » evapotranspiration efficiency « 

P 
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To a first-order approximation, we assume that the actual evapotrans- 
piration rate e^ at anytime will vary linearly with changes in the soil mois- 
ture concentrations, s. Thus, we expand J around s^ and obtain: 


®T ®T ^ 9J , . . 

™ B ~ + _ (s_a ) + . . . 

e e ds o 
P P 


From the definition of J (Equation 3.4) we can find: 


(4.5) 


C 

1 3s 


C + (2C) - - C(2C) -E • 


[t - V-t-Vj •[{-(l+Mk,)B + - B(2B)'= e] 

•e"‘^’^(l)(d+2)8^"^j' - [y [|.0e| - y f.BEjJ 

3/2r,3/2^-CE(l)8f2 - 


•e"®^E(l)(d+2)s^‘^‘‘'^J- - I -k 


^ (^ + 2) 

+ (2E) S^'' 2 ^ ^‘'(d+2) E(l) 


Then Equation (4.5) is: 


p 


(4.7) 


The following relation holds for the annual expected value of the sur- 


face runoff R ; 


-G-2o„, . -a _ 

r.a+l)a 


(4.8) 


At this point, we assume that the surface runoff occurs only during the 
storm duration. In order to obtain an expression for the average annual sur- 


face runoff rate y we can write: 
s 


where 


eTp^ ■ p 
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(4.9) 
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p ° nip /m^ nij. ^ mean storm intensity s ra, 
A r •' i 


and y = E[Rg ]/m^ 

A r 

Here again, a first-order approximation is made and the surface runoff 
function is expanded in Taylor-series around the mean 

In order to accomplish this, the derivative of the surface runoff func- 
tion, with respect to s, should be derived. In attempting to evaluate the 
derivative of the runoff function, e ^ ^^r(a+l)a a difficulty is encountered 
because the derivative of the gamma function cannot be given in closed form. 

This difficulty was overcome by approximating the function ^ = e ^ — by 

the polynomial: 

logc; = -0.806 - 1.766(loga) - 0.980(loga)^ (4.10) 

It is believed that this approximation represents satisfactorily the 
runoff function (Figure 1). 


We can now evaluate the derivative of ^ e ^ ^^r(a+l)o ^ with respect 


P 


to s , 


We find that; 




d(y„/p) 

s 

ds 


. gay 


c-1 


= e - C 


— u . 

K(l) c 


2 


+ U 


(4.11) 


where 


U = 


and 


-1.766 - 1.96(loga) 


I 67T6m I 3 '’o-' 


-4/3 


-2(1-6 )[d(l-s )l-^25-0.0375d )2-425-0.0375d^^^_^25_o.0375d) 

o o o 


^ ,1.425-0.375<1 ^ 3/3.V3 

o 


SURFACE RUNOFF FUNCTION 


FIGURE 1 


Plotted points represent: 

(logC) = -0.806 - 1.766(log0) - 0.980(logO)^ 









Tba following relation holds for the annual expected value of the ground- 


water runoff R 


ils 


E[PJ 


(4.12) 


We assume that percolation to the water table occurs during the entire 
rainy season of length m . Thus, in order to derive a representative rate 
y for the percolation, the following normalization is applied: 


«8 


E[R 1 



Ti" 


^ . 


m 


m,,m. 


o E[P^] mp 

A r 

By taking the derivative of y /p with respect to s , we obtain; 

O 


(4.13) 


^3 “ 


d(y /p) 

h 

ds 


m^mt^K(l) c s 


c-1 


(4.14) 


' s=s 




We also expand the groundwater runoff rate y around it's long-term 


average value. Thus, finally we can write to the the first approximation: 

(4.15) 


y y 

^ + C2(s-s ) 

P P 2 0 


and 


y y 

- + c (s-s ) 
p p 3 o 


(4.16) 


In order for the model to be applicable for both bare soil and for the 
presence of vegetation cover, the value of case of bare soil is 

also evaluated . 

For bare soil, the value of J (Eagleson, 1978) leads to: 
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= ~ = [- VT e"^ + [1 + ^ E] e"^ + r (|,e) ^ .|.e‘ 

, . , 2StiK(1)i|^( 1)4) (d)(d+2)s/'^^ 

- (2E)'^(eVM T-l ^ 


(4.17) 


According to the above linearizations, the conservation of water mass, 
Equation (4.1), can be written in the following finite difference form; 


(i) Rain (t £ t^) 


.7 As , 

nZr -rr = i ~ y ~ V 
At 


(4.18) 


Since surface runoff will start to be produced after time t > t^ from 
the beginning of the storm, where t^ is the time when the surface gets sat- 
urated (Equation 3.7), we must account for y in Equation (4,18) only when 

s 

> t > t . 
r o 

It also seems reasonable to assume that the percolation below the depth 


will not only be a function of the soil moisture Sj^ in the surface layer 
at time k but also of the soil moisture s^ below that layer. That is 
y = f(s, ,s ). In order to keep the equations simple, we will assume that y 


will vary with respect to the average value 


Sv + ®c 


Thus, finally we have: 


For t < t 


. f^k+1 - , r- -Pk+% 

At I ^ “ ^g S 2 


(4.19) 


For t > t 
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i - lyg P 


■[> 
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8 , -8 

k o 


<■81 ,8 Y 
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yg + SPl, 2 J 


(ii) No rain (lnter8torm period) t j< 


nZr 


f\+r\i 

r 

At j 



+ C, e (s, -8 ) - \ 

1 p k o'J \2 


yg + S » 


8, +8 
k O 


-o)] 


(4.20) 


(4.21) 


The limiting value for e„ = e,^ + C, e (s, -8 ) will be the value of the 

T 1 P K O 

potential rate of evapotranspiration e . That Impllea that e„ will be replaced 

P ^ 

by e^ until the time when the surface gets dry. 

The above equations were solved explicitey with respect to The time 

step At was taken equal to 30 minutes, i.e., we update the soil moisture every 
30 minutes. The time step was chosen to be of this order of magnitude both 
for reasons of numerical accuracy of the solutions and because this is the 
necessary time scale for conjunctive operation of the model with a GCM of the 
atmosphere. All other parameters appearing in Equations (4.19) through (4.21) 
are treated as known inputs in the model. 

Two catchments were selected to test the model, Clinton, Massachusetts 
and Santa Paula, California. They represent two contrasting climates, the 
first corresponding to a humid and the second to a se.^-.-arid region and have 
been well-studied elsewhere (Eagleson, 1978 a,b,c,d,e,f ,g) . The appropriate 
selection of parameters and necessary inputs in order to implement the model 
are presented in the next chapter. 
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CHAPTER 5 


Selection of the Appropriate Model Parameters 
5.1 Introduction 

The parameters necessary to implement the model can be divided into the 

following categories: 

a. Climatic Variables : 

HL = mean annual precipitation [cm] 

^A 


m. 


m. 


m 




= mean time between storms [days] 

= mean storm duration [days] 

= average rainy season length [days] 
s shape factor of gamma-distributed rainstorm depths 
= average annual atmospheric temperature [°C] 

= mean storm intensity [cm/day] 

JL 

= potential evapotranspiration rate [cm/day] 

The values of nu » ^ derived using the sta- 

A b r 

tistical properties of the rainfall events (Eagleson 1978b). The value 

of T is taken from measurements of the air temperature close to the sur- 
a 

face during the year. The value of e^, as it was developed in Chapter 2, 
is a function of s* t'ral climatic and surface characteristics. Here, it 
will be evaluated using a Penman-type equation. For the applications at 
Clinton, Massachusetts and Santa Paula, California, it will be set equal 
to its annual average value. In a later application, at Phoenix, Arizona 
this assumption will be relaxed and diurnal changes of will be considered. 
b. Soil parameters 

Three independent soil parameters are used in the model. These are: 
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n » effective poroalty of the soil 

2 

k(l) = saturated intrinsic permeability [cm ] 
c ° pore disconectedness index. 

By applying ecological optimality conditions (Eagleson 1982) it is 
possible, given the porosity, to estimate the appropriate values for k(l) 
and c of natural surfaces. Those conditions are described in paragraph 
(5.2). 

c . Vegetation parameters 

Vegetation is represented in the model by the percentage of vegeta- 
tion cover M and the water use coefficient k . It will again be shown 
o V 

how the value of is selected by applying ecological optimality hypoth- 
eses. Another way of obtaining is through observations, sometimes by 
using remote sensing techniques. If this is the case, then that can help 
to determine more accurately the parameter k^, as it will be shown later. 

d. Surface layer thickness 

Thn surface layer thickness is treated in the model as an independ- 

ent variable, although we know that it is a function of the root zone 
depth, and the soil and climatic characteristics of the region. Since 
the exact value of this parameter is not known and the purpose of its 
existence is to provide us with a simple conceptual model of the physical 
process which accounts for some storage of water close to the surface, it 
will be possible to fit the value of this parameter either using avail- 
able observations or solutions of more accurate numerical models. Sen- 
sitivity analysis will be performed in order to investigate the influence 

of Z on the various fluxes within the soil column, and to test the as- 
r 

sumption of many investigators, that Z must be taken equal to Ira. 

A7 


Eagleson (1982) derived several equilibrium conditions, which he hypoth- 
esized to hold in the long-term for a natural, water-limited soil-vegetation 
system. 

In a natural soil-vegetation system equilibrium stages can be considered 
to occur at different time scale. 

In the short-term it is assumed that the system tends to minimize water- 
demand stress, so, the canopy density and the plant species will take such 
values that maximize the soil moisture. That implies that the following rela- 
tions must hold for a given climate and soil: 

=0 , M = M (5.1) 

’ o 

V 



where 

s^ = average soil moisture concentration in the root zone 

M = short-term equilibrium canopy density 

o 

k = short-term equilibrium plant coefficient 
o 

Equations (5.1) and (5.2) define the "complete vegetal equilibrium". It 
was shown by Eagleson [1982], that for canopy densities > 0.42 complete 
equilibrium is not possible because Equation (5.2) cannot be satisfied. He 
further hypothesized that for a moist climate the canopy will always satisfy 
Equation (5.1) but that the species will only be in a quasi-equilibrium, so 
that the following condition is satisfied: 

k 

V 
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0 


(5.3) 



where 


E = dimensionless climate-soil parameter 

It was hypothesized that in long-term there is a synergistic symbiotic 


development of both the soil properties and the vegetation canopy which tends 
to maximize biomass production, For water limited systems, will be 

proportional to the water utilization by the plants, i.e.. 


B 'v M k e 
p o V p 


(5.4) 


For a given climate and constant k^, B^ is maximized, according to Equation 
(5.4), when is maximized. The conditions for this equilibrium then, are: 


f3M 

y 


= 0 


k(l) 


1 


= 0 


, M = M 
’ o o 


, M = M 
’ o o 


(5.5) 


(5.6) 


(8k(l) j c 

The third soil property, the effective porosity n, is assumed constant 
during this optimization procedure. 

For two contrasting climates, those of Clinton, Massachusetts and Santa 

Paula, California, and with the climatic and soil parameters given in Table 5.1 

contours of constant M for different combinations of k(l) and c were drawn. 

o 

(Figure 2 and 3) • Each contour corresponds also to a constant value of E and 

of the actual evapotanspiration . From those contours, the optimum value 

* ^ 

of the canopy density, , can be derived. That peak value of defines a 
unique pair of values of k(l) and c. Thus, under the previously developed 
hypotheses the only soil parameter that is needed to be known is the porosity 


n. 


Eagleson and Tellers [1982] provide encouraging tests of the above hy- 
potheses for different catchments. They also suggest algorithms for fitting 
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CLINTON, MA 







the values of and n, If observations of the canopy density and the total 
water yield exist. 

Also shoim in Figure 2 and 3 are curves of constant s^. It is noticed 
that the locus of optimum (maximum) s^ does not coincide with the peak value 
of The basic reason for that is that in the long-term evolution, the dom- 

inant factor is the maximization of the biomass production. Thus, although 
the system is locally optimized with respect to s^, maximization of is not 
the primary condition to be fulfilled. A more detailed discussion of this 
point is given by Eagleson [1982], 
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TABLE 5.1 


1 


Clinton, Massachusetts Santo Paula, California 


M ° 

0 

0.912 

M « 
0 

0.424 

p 

0.150 cm/day 

P 

0,274 cm/day 

" 

3 days 

S 

10,42 days 

r 

0.32 days 

t 

r 

1.43 days 

T 

365 days 

m ° 
X 

212 days 

w/e 

P 

0 

w/e “ 
P 

0 

m « 

V 

109 

% - 

15.7 


94 cm 

TOp « 

a 

54 cm 

k = 

V 

1 

k « 

V 

1 

T 

a 

8.4°C 

T 

a 

13.8“C 

< » 

0.50 

K = 

0.25 

A « 

0.578 

A 

0.0732 

k(l) = 

5.57xl0"^^*cm^ 

k(l) = 

12.27xl0"^^cm^ 

c = 

4.75 

c « 

5.25 

LThe 

values of k(l), and c 

were set equal 

to those corresponding 

to peak climatic values, according 

to the vegetal 

equilibrium hypothesis 

and the 

ecological optimality hypothesis, as they 

are described by Eagleson 


(1982).] 
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CHAPTER 6 


6 . 1 Slmulatjon of the Rainfall Process 

In order to test the model, rainfall inputs were generated, which posses- 
sed the statistical characteristics derived from historical records. Under 
the assumption of independently distributed rainfall depths, storm durations, 
and times between storms, we generated those variables with the following 
characteristics : 

Storm depth h; 

It was considered as Gamma distributed with parameters k and X. The cor- 
responding pdf was 


t„(h) 

2 2 

with = k7X, = </(X) 

Storm duration t^; 


X(X^ 


--1 -Xh 


r(k) 


( 6 . 1 ) 


It was taken as exponentially distributed with pdf of the form 

-(St 


(ty) = <Se r, ° 


( 6 . 2 ) 


where 


6 = m. 


-1 


Time between storms ty. ; 

It was also taken as exponentially distributed with pdf 


f^ (t,^) = Be“^S, t^ > 0 
% 


(6.3) 


where 


3 = m. 


-1 


5A 


The above distributions for h, t , and t, were c)iosen, because it is shown 

r b 

by Eagleson [1978b] that they can adequately represent the rainfall process. 

The generated variables, i.e., h, and t^ preserved the above-defined 
statistics. For their generation the IMSL library subroutine GGAMR was used. 
This subroutine was incorporated into the main program "Taylor. Fortran", 
which also calculates the. statistics of the generated variables in order to 
check for consistency with the historical values. (Gamma and exponentially 
distributed variables, can both be generaged with this subrountine, by making 
some slight modification of the parameters used.) 

The storm intensity was assumed uniform during the rain and thus was de- 
rived just by dividing the value of the generated h with the corresponding 
value t^. Since the storm magnitudes and storm durations were assumed inde- 
pendently distributed, the matching was performed arbitrarily by using the 
values of h and t^ in the sequence they were generated from the random number 
generator. Of course, such a matching could give rise to unrealistic values 
of i, in the extremes where independence is most invalid. However, at this 
stage, this fact will not be taken into account in testing the model, although 
its inpact should be kept in mind during the interpreted ■'.on of the results. 

The statistics of the generated rainstorm characteristics for Clinton, 
Massachusetts and Santa Paula, California are presented in Table 6.1. 

The observed differences are reasonable, since the generated variables 
which were tested, corresponded to many fewer events than the historical values 
derived from five years of observations. We also observe that the generated 
series at Clinton, Massachusetts gives an average storm depth considerably 
less than the average of five years, so we expect that to reflect in a smal- 
ler soil moisture on the average than the average annual soil moisture cor- 
corponding to five years of data. 
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TABLE 6.1 


Clinton, Mupsachusetts 



Historical (5 years) 

Generated 

(1 year) 

Stom depth 

E[h] 

= 0.86 

E[h] 

0.73 

[cm] 

Var[h] 

«=■ 1.50 

Var[h] = 

0.73 

Storm duration 

E[t^] 

■= 0.32 

?:[t ] “ 

r-" 

0.34 

[days] 

Var [t^l 

= 0.10 

Var[t^] “ 

0.13 

Time between 


= 3 

E(h] 

3.18 

[days] 

Var[t_! 

r 

=j 9 

Var[tjj] “ 

10.88 


Santa Paula, California 




Historical (5 years) 

Generated 

(1 year) 

Storm depth 

E[h] 

= 3.41 

E[h] = 

3.83 

[cm] 

Var [h] 

= 46.65 

Var]h] “ 

20.16 

Stom\ duration 

E[t^] 

= 1.43 

E[t^] 

2.34 

[days] 

Var [t^] 

= 2.04 

Var[t^] ■= 

4.69 

Time between 

E[t^] 

= 10.42 

E[tj^] = 

11.90 

[days] 

Var [t^l 

= 108.57 

Var[t^] = 

66.30 
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CHAPTER 7 


Presentation of Results 

7 . 1 The Evapotransplration, Surface Runoff, and Percolation Functions 

Before applying the model described in Chapter 4 for simulating the soil- 

moisture concei.tratlon during the rainy season, it is essential to present the 

rate functions of evapotranspiration e , surface runoff y , and percolation y , 

T ^ 8 

which will be linearized around the average annual soil moisture, s^. 

The actual evapotranspiration rate e^ is given by Equation(3.6) , where it 
appears through the expression of evapotranspiration efficiency J, i.e., nor- 
malized with the value of the average annual (seasonal) potential evaporation 

rate e . By using the climatic variables and the soil and vegetation paramet- 
P 

ers derived under climatic climax conditions and shown in Table 5.1, for the 
catchments of Clinton, Massachusetts and Santa Paula, California, the evapo- 
transpiration efficiency J can be plotted as a function of the soil-moisture, 

s. 

Figure 4 shows the J(s) function for the bare soil case (M^ = 0) and Fig- 
ure 5 shows J(s) when the presence of vegetation is taken into account. Also, 
shown in Figure 4 is the evaporation efficiency function that is used by 
Manabe [1969] in his GCM. This follox-js a linear Budyko-Type parameterization 

for which J = 1 if s > s, and J = — if s < s, , where s, is a critical value 

Ic ^ ^ ^ 

of s defined by Sj^ = 0. 75 x and is the degree of soil saturation with- 

c c 

in a soil layer from the surface to Im depth, corresponding to the field capa- 
city 0^ . The value of the field capacity used in Manabe ’s General Circulation 
c 

Model is held uniform over all areas of the Earth and is set equal to 15 cm. 
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The difference in the J(s) functions for Clinton and Santa Paula can be ex“ 

plained by the different climatic and soil properties of the two catchments* 

Clearly, Manabe’s J(s) function overestimates the evapo transpiration in both 

cases. In Figure 5 we observe the apparent influence of the vegetation cover 

M on the evaporation efficiency. In this case according to the assumptions 

developed in Chapter 3, the limiting value of e„ becomes M k e . For the hu- 

T o V p 

mid climate of Clinton, we expect that the deviations of s around the average 
annual value s^, will not be very high and thus we will always operate in the 
region where J = 1. Otherwise, the linearization procedure described in Chap- 
ter 4 will not give accurate results. With a value of s = 0.72, as derived 

o 

from the annual water balance, the above assumption is very reasonable for the 
humid climate. For the semi-arid climate of Santa Paula, we observe that the 
function J(s) is very close to a linear form, for values of s in the neighbor- 
hood of - 0.55. That implies, that the use of a linearized function of J 
around s^, will give accurate results for this case. It must be pointed out 
that when a constant value of e^ is used, equal to its annual average value e^, 

then the actual evapo transpiration rate e_ will tend to e as s increases above 

T p 

the value of Thus, this function is expected to give fairly accurate re- 

sults when applied in a real case of successive rainy and dry periods, if e^ 
is supposed to be held constant. But if e^ is changing, as would be the case 
in reality, then the value of e^ obtained from this function will tend to the 

value of e whenever the surface becomes saturated and not to the actual value 
P 

of e . Thus, if a changing value of e is to be used, the time t from the 

beginning of the interstorm period until the surface gets dry must first be 

calculated. Before that time, evaporation will occur at the actual potential 

rate e and after that time control will pass to the soil and the value of e 
^ 60 


obtained through the previously described linearization procedure will be used. 

The average annual surface runoff rate y is given by; 

s 

-G-2o 

y^ « e r(a+l)a“°iOp (7.1) 

r 

In Figure 6 y is plotted versus s, for the catchments of Clinton and 
s 

Santa Paula. 

The average annual percolation rate y to the water table is given by; 

S 

yg = KCDs*" (7.2) 

The function y versus s is shown in Figure 7 for both Clinton and Santa 
Paula. 

The functions y (s) and y (s) indicate high non-linearities between those 
3 g 

fluxes and the soil moisture s. Thus, we must expect that a linearization 
around the corresponding average soil moisture s^, for each climate, c-^n intro- 
duce errors in the estimation of those rates, especially when the deviations 
from the average value become high. This fact should be considered with at- 
tention to the interpretation of the results. That is, since those functions 
appear in this case to be convex, we expect to underestimate the surface runoff 
and percolation rates, whenever s > and overestimate them when s < s^. 
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Surface Runoff Function 
FIGURE 6 
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(cm/day) 




7.2 Simulation of Soil-Moisture Concentration During the Rainy Season Using 


a Constant Value of ep ° Cp 

Equations (4.1A) and (4.16) were solved for each time step (At = 30 
minutes), with respect to using the generated rain storm events. The 

depth of the surface layer which accounts for storage was treated as an in- 
dependent variable, i.e., many simulation runs were performed with different 
values of in order to observe the sensitivity of the fluxes with respect to 
that parameter. The capillary rise from the water table was taken equal to 
zero, assuming that the water table was deep enough that it did not have any 
impact on the fluxes occuring close to the surface. The climate and soil pro- 
perties for Clinton and Santa Paula were those presented in Table 5.1 which 
were derived using ecological optimality hypotheses (Eagleson 1982). 

The computer program named "TAYLOR. FORTRAN" was set up to perform a simu- 
lation of the soil-moisture concentration in the surface layer. A complete 
description of this program is given in Appendix 2. 

The soil-moisture concentration, s, as a function of time, for Clinton, 

Massachusetts is shown in Figure 8. Two different cases are presented; one 

with Z = 100 cm and one with Z = 50 cm. As is expected for the case where Z 
r r r 

is smaller, the soil moisture fluctuates over a larger range. The results sho\im 
in Figure 9 correspond to a vegetation cover M^ = 0.912, which is the climatic 
climax value (Eagleson and Tellers, 1982). When bare soil was assumed (M = 0), 
there was no change in the results, because in the humid climate of Clinton, 
evaporation from the bare soil occurs at the potential level (climate control- 
led) and because we have taken k =1, its optimum value (Eagleson and Tellers, 

V 

1982). 
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Analogous results of s versus time for Santa Paula, California and • 

100 cm and ° 50 cm, are shown in Figure 9. We observe that at Santa Paula 
we have larger deviations of s around the mean s^ compared to those at Clinton. 
This is due to the much longer interstorm periods and longer storm durations 
of the climate of Santa I aula. The results shown in Figure 10 correspond to 
a vegetation cover ° 0.42A (optimum value). When M was set equal to zero 
(bare soil), small differences occured in the soil-moisture concentration. 

This was due to the fact that here also the evaporation from the bare soil is, 
on the average, rather high (J(s^) = 0.84) and again “ 1. 

A sensitivity analysis was performed for the humid climate of Clinton and 
the semi-arid climate of Santa Paula, with respect to the value of the parame- 
ter Z^. The results are shown in Figure 10. The horizontal axis corresponds 
to values of ranging from 20 cm to 120 cm. On the vertical axis, the cumu- 
lative yield and cumulative evaporation at the end of the rainy season are 
plotted. For Santa Paula, California, it is interesting to observe that there 
is a range of values of Z^ from approximately 60-120 cm, where those fluxes 
are insensitive to the value of Z^. When Z^ becomes small enough, evaporation 
is rapidly reduced. This is due to the fact that soil moisture is exhausted 
very fast during the interstorra period, the surface becomes dry faster and con- 
trol passes to the soil earlier than before. On the contrary, yield increases 
rapidly with small values of Z^, because in that case, during the rain, the 
soil-moisture concentration in the surface layer increases very fast and the 
surface becomes saturated at earlier times. Thus, surface runoff is produced 
more frequently aud at earlier times during the rain. 
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Simulation of Soil Moisture During the Rainy Season 
(Santa Paula, California) 

FIGURE 9 
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Sensitivity of Fluxes to Surface Layer Thicl}.ness 


FIGURE 10 
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CUMULATIVE EVAPORATION AT THE END OF THE RAINY SEASON (cm) 



For Clinton, Massachusetts, the cumulative yield docs not seem to become 
very sensitive to the value ol 2^, when the latter becomes small. For values 
o£ ranging from 0.4-J.m, there is a maximum change of 3cm in the cumulative 
yield. The rote of rhang'' becomes much smaller when exceeds Im. Cumulative 
evaporation was fuune to remain constant for all values of Z^, indicating that 
the exfiltration process was always under climate control. At this point, it 
must be noted that, because of the model's structure, this sensitivity analy- 
sis becomes Invalid for the humid climate of Clinton when is low and the 
value of s is such that control must pass to the soil. That is, the lineariza- 
tion around the value of the average annual soil Eolsture s^, is not valid in 
this case of e^ << e^. However, for a humid climate, we can say apriori that 
should not be very low and so assume that soil control does not occur. Othe 
wise, it would be necessary to change the linearization procedure, and thus 
reduce the general applicability of the model. 

It must be noted that in all cases examined for both climates, the mois- 
ture in the surface layer was never completely exhausted. would be ex- 
hausted however, if even smaller values of Z are assumed. However, such Z 
are physically unrealistic because Z^ is defined to include all exchangeable 
moisture. 

We also observe that, because of earlier passage to soil-control when 
M 5 * 0, the cumulative evaporation is greater when M 0, for most values of 
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7.3 Comparison With A Numerical Model 


In order to test the predictive capability of the model developed In this 
study, the soil-moisture concentration and the water fluxes obtained from it, 
were compared to those obtained by an "exact” numerical model. The numerical 
model used was that developed by Milly and Eagleson [1980]. This model asaximes 
a one-dimensional representation of the physical system and solves the coupled s 
non-linear partial differential equations governing mass and heat transport in 
the soil, using the Galerkin finite element method. In the present comparison^ 
an isothermal version of the model was applied. The vertical soil column was 
taken equal to 5m and the influence of the water table was considered to be 
negligible. A constant flux (K(0^) = constant) boundary condition was assumed 
at the bottom. The surface boundary condition was changed according to the sur 
face moisture state. During precipitation* infiltration takes place at the 
precipitation rate, until the soil surface reaches saturation. After that hap- 
pens, ponding of water on the surface occurs, thus prorlucing surface r^ .-off . 
During evaporation, the evaporation rate is equal to the potential value until 
the surface becomes completely dry and control passes to the soil. After this 
time, evaporation proceeds according to the exfiltration capacity of the soil. 
The values for the precipitation intensities, storm durations, and interstorm 
periods used were those obtained from the aimnlated rainstorm events, as de- 
scribed in Chapter 6. The value of the potential evapottanspiration was as- 
sumed constant throughout the simulation period and equal to it’s annual aver- 
age value, and thus, as discussed in Section 7.1 it wac? not necessary to cal- 
culate a time t^ during the evaporation period. 

The computer program SPLASH .FORTRAN developed by Milly [1980] was used in 

order to obtain the numerical solution. Many runs were performed, varying 
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the number of finite elements within the soil column and several other para- 
meters, in order to achieve convergence of results. This procedure is de- 
scribed with more details in Appendix 2, The results from this comparison are 
shown in Figure 11-16, and correspond to solutions of the numerical model where 
convergence was achieved. Figures 11-13 show the storage change, the total 
yield (surface runoff and percolation) and evaporation flux respectively, for 
Santa Paula, California and for 10 consecutive simulation periods. Each simu- 
lation period corresponds to either a storm or an interstorm period, with in- 
tensities and durations as generated by the procedure described in Chapter 6. 

The storage change represents the deviation from the initial soil-moisture con- 
centration. From the results of the numerical model, it was found that the sur- 
face became dry at the eighth simulation period v^hich implies that control pas- 
sed to the soil at that time. 

It is observed that the differences between the two models are not big 
and never exceed 1 cm storage. The analytical model follows the numerical 
solution very faithfully. The water fluxes outside the surface layer are shown 
in Figure 12. Again, the differences between the two models are very small. 

The evaporation flux is shovm in Figure 13. The differences here do not ex- 
ceed 0.5 cm. 

Figure 14-16 show the same comparison for Clinton, Massachusetts, and for 
20 simulation periods, generated as discussed before. From Figure 14, we ob- 
serve that the storage change is consistently -^er estimated by the analytical 

model, but again the differences between the tv* relatively small and never 

exceed 0.6 cm of storage. In Figure 15, the total yield at each simulation 
period is plotted and the agreement between the two models is considered as very 

satisfactory. Evaporation fluxes are shown in Figure 16. For the humid cximate 
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ANTA PAULA, C A 
YTICAL MODEL (Zr = I00cm) 
EPICAL MODEL (41 nodes) 
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SIMULATION PERIOD 

Ion Produced by the Analytical and Numerical Model 



of Clinton, evaporation was always at the potential level and thus there is no 
difference between the two models. Thus, all differences in storage change 
can be explained by the differences occuring in the prediction of the yield. 


In general, it should be noted that the solution obtained by the analy- 
tical model is in very satisfactory agreement with the numerical solution 
for both climates. Of course, it must be kept in mind that the tested ver- 
sion of the numerical model was a simplified one, since Isothermal conditions 
were assumed. 

7.4 Comparison with Manabe's [1969] Parameterization 

Manabe's [1969] landsurface parameterization was also compared with the 
model developed in this study. 

Manabe [1969] uses the concept of field capacity s^ in his soil-moisture 

c 

model. He assumes a surface layer of Im and defines a critical value of soil- 

moisture Sj^ given by: « 0.75 x . Then he assumes the following equa- 

c 

tions to hold for the water and vapor fluxes at the surface: 
i. Evaporation 

if 3 > 5^, 


if s < s^, e^ 


s 

e . — 
P 


where 


e^ is estimated from an aerodynamic type equation 
ii. Soil moisture 

3s 

if s = s- and i > e , ir- =» 0 and y == i - e 

p ■'s p 

and if s < = i - e^ 

The value of the field capacity chosen by Manabe for all applications 
was 15cm, which for a soil layer of Im and a porosity of 0,35 corresponds to 
a soil-moisture concentration given by = 0.A2. 


78 


Using the same initial condition (0 “ 0^ ) for both models, the stor- 

c 

age change versus the simulation period, is shown in Figures 17 and 18 for 
Santa Paula, California and Clinton, Massachusetts. It is observed that 
there is a big difference between the two, which is on the order of 4cm for 
Santa Paula and 3cm for Clinton. Since the first simulation period corres- 
ponds to a storm and the initial condition of soil moisture is set equal to 
the field capacity, it is expected that Manabe's model will not produce a 
storage change during that period because according to his assumption, soil- 
moisture cannot exceed the value of field capacity. If the first simula- 
tion period was an interstorm period, the storage changes produced by both 
models would not have so much difference. But big differences in storage 
will occur later on, when a precipitation event comes and soil-moisture 
reaches the value of the field capacity. 

Since good agreement between the presently developed analytical model 
and the "exact" numerical model has already been established, it appears 
that Manabe's model does not represent the system very well. His model fails 
to capture the time variations of yield. This is shown for Clinton, Massa- 
chusetts in Figure 19. Values of the yield for Santa Paula obtained using 
Manabe's model are not shown graphically here. It was found, however, that 
a total yield of 8.7 cm was predicted for Santa Paula during the first sim- 
ulation period by this model, while the yield was zero for all other simu- 
lation periods. As is cHown in Figure 12, this is much different from the 
value predicted by both the numerical model and the analytical model devel- 
oped here. Better results could possibly be obtained if a different value 
for the field capacity was chosen for Manabe's model, but what means can be 

used to evaluate this field capacity a priori? 
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SIMULATION PERIOD 

Comparison of Storago Change Produced br Manabe's Model and the Analytical Model 
(Clinton, Massachusetts) 

FIGURE 18 
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In general, the analyti^a:. model suggested in this study seems to be an 
improvement to the landsurface parameterization, mainly because it is simple, 
physically based and gives consistent and accurate results when compared to 
numerical solutions. 

7 • 5 Soil Moisture Simul .- :n with Changing Value of ep 

The model was also tested with real measurements of soil-moisture concen- 
trations obtained from an experimental field at Phoenix Airport, Arizona. 

Values of the meteorological variables were available every half-hour, so that 

is was possible to estimate a changing value of e , using either a Penman-type 

P 

equation or an aerodynamic equation. More precisely, the following measured 

aata were available: Net radiation at the surface, Air Temperature T^, 

wind speed U , and vapor pressure e at screen height, surface temperature 
a o 

T at depth of 1 cm and average soil-moisture concentration in three layers 

O 

below the surface (from 0 - 10cm, from 10 - 50cm, and from oO - 100cm). The 
data corresponded to seven days of measurements from 5 March - 11 March, 1971. 
Details of the experimental field and measurement procedures are given by 
Jackson [1976]. Briefly, the soil consisted of Adelanto Loam, was reasonably 
uniform to about 100cm and had been cultivated numerous times during past years. 
The soil properties for the Adelanto Loam and the climatic variables at 
Phoenix Airport, are given on Table 7.1. A graph of hydraulic conductivity and 
diffusivity as a function of soil moisture 0 is given by Jackson [1976]. Be- 
fore the experiment took place, the field was irrigated and during the seven- 
day period of measurements, no precipitation was measured by rain gages at the 
Phoenix Airport, although some traces did occur. 

In order to compare the results of the model with the measured values of 

soil-moisture, the latter were averaged over the Im surface layer depth. 
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TABLE 7.1 


Phoenix Airport - Arizona 

6p « 0.263 cm/day 

m = 7.27 days 

ip.j. = 0.11 days 

r 

= 365 days 

K = 0.50 

w/e - 0 

P 

T =21.3 °C 

a 

m =45 

V 

iiL = 19. '5 cm 

A 

n = 0.35 (assumed) 

k(l) = 2.68x10 ^ cm^ 

c = 6.5 

[The climatic variables shovm in this table were derived by using the com- 
puter programs HODCOP and RAINSTAT developed by Restrepo and Eagleson (1978) 
for the interpretation and analysis of NOAA hourly precipitation data tapes. 
The soil properties for the Adelanto loam are given by Jackson (1965, 1S79) 
and Mualem (1976).] 
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The value of the potential evaporation rate e^ was first estimated from 
a Penman-type equation, as described in Chapter 2, where the surface heat flux 
was neglected and the surface roughness was taken equal to 0.03cm. 

As is also mentioned by Jackson [1976], the surface became dry after the 
fourth day of the experiment. This implies thcii; after that time, flux control 
passed to the soil. Since the model used accounts only for soil moisture with- 
in the bulk volume and since it uses an evaporation efficiency function de- 
rived using the value of the annual (or seasonal) average evapotransplration 
rate e^, it cannot accurately locate this change, especially when e^ is much 
higher than e^, as it was also discussed in Section 7.1. To surmount this 
problem, it was necessary to calculate apriori the average time t^, until the 

surface becomes dry. This time is given by; (Eagleson, 1978d) 

„2 


" 2S 2 


(7.3) 


where S is the exfiltration "desorptivity" defined for a dry surface by; 


S = 2s 
e o 




m 


(7.4) 


and M was assumed equal to zero. 

Using the values of t}ie parameters as defined in Table 7.1, it was found 
that t^ = 3.92 days, which is very close to the value of four days mentioned 
by Jackson. 

Equation (4.16) was now solved as before. The results for the soil- 
moisture concentration in the layer of Im, are shown in Figures 20 and 21, 
evaluated using two different values of e (the annual average value e = 0.263 

P Pa 

cm/day, and the average value of during the seven-day period, equal to 0.323 
cm/day) . It is observed that the results are extremely encouraging and also. 
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FIGURE 20 
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as expected, the fitting is better when e^ is the actual average for the seven- 
day period. 

Some of the small discrepancies in these comparisons may be due to the un- 
gaged (trace) precipitation that occured during this period. 

In order to perform the simulation described above, the computer program 
ARIZ. FORTRAN was prepared. It uses as an input file the available meteorologic 
and soil-moisture characteristics, obtained every half-hour. 

It must be noted that the ecological optimality hypothesGs were not ap- 
plied for soil-parameterization in the Arizona experiment for two reasons. 
First, because the field has b ,en cultivated and is thus not in its natural 
state, and second, because soil properties were available from measurements. 

Evaluation of the soil moisture concentration and the evaporation flmces 
during the seven-day period will be shown in the following Section 7.6, vrhere 
the thermodynamic coupling will be applied. 

7 . 6 The Thermodynamic Coupling 

The developed soil-moisture model was operated conjunctively with a ther- 
mal balance model, in order to estimate the surface temperature T^. Two meth- 
ods were tested, using the Arizona data. One used the force-restore method 
[Deardorff, 1978] and another used the thermodynamic equilibrium equation 
[Edlefsen and Anderson, 1963]. 

a. The force-restore method. 

The force-restore method was described in Chapter 2, but the basic equa- 
tions used are repeated here for convenience. Written in finite difference 
form, as they were solved in the present study, those equations take the form; 
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k+1 .J1 

J5 . 2tt k ^ k 

PsS^l '^l g 2 ^ 


(7.5) 


G = R - H - LE 
n s 


(7.6) 


rnk+1 - k 
2 ^ 2 


= G/(p^c^d^) 


(7.7) 


- ground temperature at the surface (®K) 

2 

= heat flux into the soil [cal/cm .sec] 

= net radiation at the surface [ly/sec] 

2 

= sensible heat flux [cal/cm .sec] 

9 

= water vapor flux [g/cm"sec] 


= 86400 sec 


T 2 = mean soil temperature over layer of depth d 2 (°K) 

L = latent heat of vaporization [cal/g] 

c^ = specific heat of soil [cal/g. °K] 

3 

p == density of soil [g/cm ] 
s 

d^ = soil depth influenced by the diurnal temperature cycle [cm] 

d 2 = soil depth influenced by the annual temperature cycle [cm] 

l< 

The value of d is given by: d- = (k t,) , where k is the soil ther- 

X X S X s 

mal diffusivity. Assuming a soil porosity n = 0.35, the volumetric heat capa- 
city of the soil pc is estimated by (DeVries, ”Heat Transfer in Soils"), 
s s 

P_c = (1-n) 2x10^ +0 (4.2) x 10^ + (n-0)c 

s s a 

where 


G = volumetric water content 


heat capacity of the air 
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Using the Arizona data, we obtain: 

pc = (1 - 0.35)2x10^ + 0.0805x4.2x10^ + 0.269 x 1.25 x 10^ = 1.6384369 

s s 

X 10^ 


or p c = 0.39 cal/cm^°K 

s s 

The thermal conductivity of the soil X is obtained using the value given 

by deVries for loam and for n = 0.35. It is found that X = 1.44 cal. cm ^sec “C 

Thus, the soil thermal diffusivity will be given by: 

k = — ~ — = 0.013 cm^/sec. 
s pc 
s s 

and 

d^ = (0.0x3 X 86400)^ = 33.51 cm. 

This value is close to the value of 31,89 applied by Lin [1980] for the same 
experimental field. 

The value of is given by; 
d„ = (365 k T,)^ = 640.21cm 

The sensible heat flux H was evaluated from the equation suggested by 

s 

Anderson [1976]: 


H = p c ,C„,U (T - T ) 
s a p U a g a 


cal 


sec-J 


(7-8) 


Under climate-controlled conditions where E = E , the water vapor flux 

P 

was evaluated by the aerodynamic relation: 


p x0,622 


E = 
P 


-.C U (e - e ) 
was a 


^ 


2 

cm sec 


(7.9) 


In these equations: 
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e = saturation vapor pressure at surface temperature 
s 

e = vapor pressure of the air at screen height 
a 

= wind speed at screen height 

T “ air temperature (°K) at screen height 
a 

= density of air 

Cp ■ specific heat of dry air 
“ atmospheric pressure 

C„ and C are coefficients equivalent to the drag coefficient for sen- 
sible heat and water vapor flux respectively. Under neutral conditions those 
are given by (Anderson 1976): 






where 




(7.10) 


k = Von Kerman's constant (0.4) 

Z " screen height 

a 

= surface roughness 

Deardorff [1968] computed the ratio of each of those coefficients to its 
value under neutral conditions, and his results are described by: 




w 




1.0- (An^)"^ 


+ i 
^ 2 


.-n“l 


„ , 1+X V , o ft / 1+X^ 

Jln(— 2 ~ ) +2in(-y-) 
2 - 1-1 


2tan ^ (x) 


. fl - 2Jln(2^)'- 


(7.11) 


where 


(1 - 16 ^) 


1/4 


and L is the Monin-Obukhov length which can be related to the bulk Richardson 

number (R, )_ through the relation: 

1 B 
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(7.X2) 


Z 

a 

L 


U.Oh/(C„)^ • <Vb 


Z 2 

[1 - [an(i^) + 2to(^) - 2tan"^(x) + f]"^ 


where 


‘Vb 


2S-ZJVT^) 
«g«a> '“a" 


(7.13) 


Using equations (7,10), (7.11), and (7.12), a table of corresponding values 
of the ratio of bulk transfer coefficients C„ and C to their neutral value 
for different Richardson numbers (R.)„ and for different values of surface 

1 D 

roughness Z^, can be constructed. Such a table is shown by Anderson [1976, 

page 19]. That kind of table was used intiie analysis performed here in order 

to determine the transfer coefficients to be used. 

Equations (7.5), (7.6), and (7.7) were solved simultaneously with the 

soil-moisture Equation (4.16). At each time-step, which was equal to 30 min- 

k+1 

utes, new values of for the surface temperature and of for the soil- 

moisture, were estimated. The parameters of the surface roughness Z^ and of 
the initial deep soil temperature T 2 were varied in order to obtain the best 
fit with the measurements of surface temperature and soil moisture. 

The changing value of e^ was evaluated through the use of the aerodynamic 
Equation (7.9). This value is used until control passes to soil. 

The results are shown in Figures 22 through 29. In Figure 22 the surface 
temperature is plotted and compared with the solid line which represents the 
measured values. The transfer coefficient was set equal to (Cjj)j^ *= 0.00277 
and the initial deep soil temperature was set equal to T 2 “ 11°C. The fit- 
ting can be considered as satisfactory, although we observe that for the firs'. 
60 hours the surface temperature is overestimated by the model at the peaks 


and after that point it is underestimated at the peaks. 
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In Figure 26, the daily evaporation rate obtained from the model is com- 
pared with that measured by lyaimeter (Jackson, 1976). It is seen that the 
calculated evaporation from the first day is less than the measured one. That 
could explain the overprediction of surface temperature observed in Figure 22, 
at the first day. That is, the higher actual evaporation makes the surface 
cooler than that predicted from the model. 

Another fact that must be mentioned is that the measurements of surface 
temperature are at a depth of 1cm below the surface. Since the model assumed 
evaluates the temperature exactly at the surface, a discrepancy between the 
two could be justified. J. D. Lin (1980) mentions that high temperature gra- 
dients, as high as a difference of 20‘’C in 2cm, can occur near the ground sur- 
face during most of the daytime, which supports what was said before. 

Figures 24 and 25 show the results obtained using the same value for the 
transfer coefficient (C„)., = 0.00277 as before, but with a different initial 
value for the deep temperature T„ = lA^C* We observe that soil-moisture con- 
centration is not predicted as well as before. 

It has been argued by Bhumralkar (Deardorff, 1978) that T 2 can be esti- 
mated as the average air temperature during the previous 24 hours. If this 
argument is correct, it is possible that an initial value of T^ = ll^C, al- 
though it seems low for the Phoenix climate (where the annual average air tem- 
perature is about 21 “O could indeed have occured. 

Figures 27 and 28 show the results of the comparison, when a larger 

transfer coefficient (C„)„ =■ 0.0057 is assumed. Clearly, for this high value 

H N 

of the soil-moisture concentration (Figure 28) is very much underpre- 

dicted by the model. 
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A complete sensitivity analysis of the errors in the predictions of soil 
moisture and surface temperature with the values of surface roughness and in- 
itial deep soil temperature is shown in Table 7.2. It seems that a value of 

Z « 0.05cm and of = 11°C gives us results which predict fairly accurate- 

i 

ly both surface temperature and soil moisture. This can be confirmed both 

from Table 7.2 and from Figures 22 and 23. From Table 7.2, it is observed 

that combinations of and T^ with even small^.r errors do exist, but the 

differences are very small compared with the errors obtained when Z^ = 0.05cm 

and ~ 11°C are selected. The problem of a priori selection of Z^ and 

T^ remains however, 
i 

It should be specially noted that when control passes to the soil (after 
the fourth day), the prediction of surface temperature is very accurate, which 
implies that the analytical model developed here for soil moisture fluxes, 
gives reasonable estimates of the actual rate of evaporation. 

Manabe*s model cannot be compared to the analytical model developed here 
for the Pheonix, Arizona experiment because the value of the soil field cap- 
acity, , is not known. 

b. The thermodynamic equilibrium equation. 

The thermodynamic equilibrium equation is given by: 

e ^g(s,T )" 

TW) = “rt 

s g L g J 

where e is the vapor pressure at the soil surface. 

Equation (7.1A) was developed by Edlefsen and Anderson (1943) and im- 
plies that the water and vapor are in thermodynamic equilibrium. It has the 
attraction of involving only known variables and thus does not require esti- 
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TABLE 7.2 


Fhoenlx, Arizona 





Cumulative Error 



Cumulative Error 

of surface 

Surface 

roughness 

Z (cm) 

0 

Inititai deep 
soil temper- 
ature T„ (°C) 

of Soil moisture 

nredictions 

MES CAL| 

^ '®i “ ®i 

i=0‘ ' 

Temperature 
predictions 
^37 MES CAL 

i=0' ®i ®i 

0.5 

13 

4.898 

563.78 


11 

3.815 

629.19 


9 

2.738 

683.88 


7 

1.935 

1008.22 


5 

1.063 

1327.97 


3 

0.457 

1696.69 

0.05 

15 

1.300 

724.60 


13 

0.541 

610.21 


11 

0.582 

571.49 


9 

1.373 

628.97 


7 

2.155 

824.57 


5 

2.790 

1135.92 

0.1 

15 

2.339 

639.97 


13 

1.422 

561.12 


11 

0.535 

550.40 


9 

0.547 

627.94 


7 

1.381 

842.12 


5 

2.016 

1162.72 

0.25 

15 

4. 1.10 

539. -4? 


13 

3.250 

519.76 


11 

2.253 

565.17 


9 

1.278 

685.42 


7 

0.474 

920.04 


5 

0.550 

1254.45 
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mates such as T 2 and of the force-restore method. In order to apply it 

during the exfiltration process, it must be assumed that a quasi-steady state 
thermodynamic condition is reached at each time-step, when e^ is calculated. 

It must also be noted that Equation (7.14) ignores the influence of the ad- 
sorptive force-field, which can become important for a dry soil. In order 
to apply Equation (7.14), the dependence of i|; on T was ignored, assuming 

o 

that the Influence of T on li/ is not important and that the primal variabil- 

S 

ity of ij; comes from variations in soil moisture. 

Also (7,14) was applied only for the case where the surface is dry. 

When the surface is wet, the surfac»i termperature was estimated by using again 
the force-restore method. If instead of doing so, it was set equal to the 
air temperature, a big discrepancy between the measurements and predictions 
would have been observed. 

At each time step a value of the actual evaporation was determined by 
the model. By using the aerodynamic equation, a value for the vapor pressure 
e at the surface was calculated. From the soil-moisture model, the value of 

ip was estimated by using the current value of s. The value of the surface 

Ic h 

temperature T of the previous time step k was used to evaluate (RT ) since 
8 S 

this factor is not very ^snsitive to changes at T , Then, Equation (7.14) 

was solved for since the only unknown now was e^ (T^"*~ ), which is a 

g ’ s g " ’ 

k+1 

function of T 

g 

The results are shown in Figures 29 and 30. It does not seem that the 
surface termperature is well estimated compared to the results of the force- 
restore method. After the fourth day when the surface dries and the thermo- 
dynamic equation starts to apply, the surface temperature is systematically 
underpredicted. 
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CHAPTER 8 


Summary, Conclusions, and Reconunendatlons for Further Research 


8.1 Summary 

In this study^ a simple analytical model was formulated to parameterize 
the water fluxes at the landcurface. A short-term water balance equation was 
solved during precipitation and interstorm periods, assuming all exchange of 
moisture to take place within a single surface layer. The evaporation and 
yield fluxes were assumed to vary linearly around their annual average values, 
as given by Eagleson (1978). 

Successive rainstorm events and interstorm periods were generated in 
order to test the model. Soil-moisture concentration, within the surface layer 
was predicted every half-hour. The storage change and the evaporation and 
yield fluxes obtained from the analytical model during long simulation periods, 
were compared with those obtained from a numerical model (Milly 1980) and with 
a simple parameterization model (Manabe 1969) currently used in GCM's. This 
was done for two contrasting climates, chose of Clinton, Massachusetts and 
Santa Paula, California. 

Finally, the analytical model developed here for soil-moisture fluxes was 
operated conjunctively with thermal balance models, in order to predict the 
surface temperature. Results of the obtained soil-moisture concentration and 
surface temperature for this latter case, were compared with available measure- 
ments. 

Two cases, one in which the potential evapotranspiration rate e^ was held 
constant at its annual average value and one in which e^ was allowed to change 
with time, were considered aad the necessary modifications of the model for 
each case were discussed. 
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Sensitivity analysis was perfornied with respect to the depth of exchange- 
able moisture within the soil and also with respect to the surface roughness 
and deep soil temperature, when the thermodynamic coupling was considered. 

For the catchments of Clinton, Massachusetts and Santa Paula, California, 
the soil and vegetation properties used were those obtained from the appli- 
cation of ecological optimality hypotheses (Eagleson 1982) , According to these 
hypotheses only effective soil porosity n and the climate are necessary in 

order to determine k(l), c and the optimum vegetation properties and . 

^0 

In order to apply the analytical model developed here to determine the 
landsurface boundary condition for use in GCM’s of the atmosphere, the follow- 
ing steps should be followed: 

1. Obtain at each grid point on the landsurface the representative cli- 
matic and soil parameters necessary to implement the model. Those para- 
meters are described in Chapter 4 and the use of ecological optimality 
hypotheses, In order to reduce their number is also discussed. 

2. Make an estimate of the surface layer thickness. 

3. Estimate the average storm intensity and average potential evapora- 
tion rate every half-hour (or appropriate At) according to whether it is 
a precipitation or an interstorm period, respectively. If there is a 
precipitation period, calculate the time t^ until the surface becomes 
saturated and let surface runoff be produced after that time if rain con- 
tinues. If there is an interstorm period, calculate the time t^ until 
the surface becomes dry and calculate evaporation before that time by 
setting it equal to the value of the (changing) potential evaporation 
rate . 
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A. SolVtf the linearized Equations (4.14) and (4.15) if a storm has 
occured, or Equation (4.16) if it is on interstorm period, every time 
period. Thus, an updated value of soil-moisture concentration will be 
obtained. Updated values of actual evaporation rate and of yield rate 
will also be obtained. 

5. If the surface temperature is to be calculated, Equations (4.14)- 
(4.16) can be solved conjunctively with the equations of the force-restore 
method. In order to do that, estimates of the surface roughness and the 
initial deep soil temperature are necessary. In addition, knowledge of 
several meteorological variables at each time step will be necessary in 
order to estimate the changing value of e^, which influences the ther- 
mal and water balance equations at the surface. 

8.2 Conclusions 

The model was tested using simulated rainstorm events for two conti.'JSting 
climates and in both cases it was found to agree reasonably well with the sol- 
ution of the numerical model. It was also tested against real measurements 
of soil-moisture during an evaporation period, and again it was found that it 
made very accurate predictions. 

Thus, from the results obtained in this research, it can be stated that a 
simple second-order Budyko-type parameterization of the landsurface, compares 
favorably with "exact" numerical solutions for exchanges of water through the 
surface. Also, the parameterization suggested here is an improvement over 
the first-order Budyko-type model of Manabe (1969) . 

A range of depths of the soil-moisture layer was found for both tested 
climates, within which the cumulative evaporation and yield fluxes were rather 
insensitive to the layer depth. This range appears to include the actual 


root-zone depth. 
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In order to estimate the surface temperature, it was found that the force- 
restore method is superior to the application of the thermodynamic equilibrium 
equation, but again difference with real measurements, although small, u’d 
exist. It must also be pointed out that in order to apply the force-restore 
method conjunctively with the analytical soil-moisture model developed here, 
parameters such as the surface roughness and the initial deep soil temperature 
must be either known or fitted apriori to available data. 

It should also be said that the very good agreement between analytical 
and numerical solutions for Santa Paula and for Clinton was obtained using 
soil properties derived from ecological optimality hypotheses. This provides 
one more indication of the applicability of these hypotheses. 

8 . 3 Suggestions for Further Research 

Using as a basis the simple landsurface parameterization developed in 
i Uxs study, the following additional studies should be carried out; 

1. Test the model at other catchments particularly under soil-controlled 
conditions and for longer simulation periods. 

2. Compare the model with other simple pararaeterizations in addition to 
the one suggested by Manabe (1969) . 

3. Investigate cases where 5 ^ 1 and possibly test the analytical mod- 
el with more accurate niomerical models which include vegetation. 

4. Investigate the relation between the soil-moisture layer thickness 
and the soil properties k(l) and c for different climates. 

5. Further investigate the sensitivity of the force-restore method with 
respect to the surface roughness coefficient and the deep soil tempera- 
ture.. 

6. Use the short-term water balance model developed in this study con- 
junctively with measurements of soil moisture (obtained for example by 
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remote sensing). Thus, one equation representing system dynamics and one vec 
tor of observations will bo available to apply optimal linear estimation te^h 
niquoB (linear filtering) and to make predictions of soil moisture. 


no 
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APPENDIX 1 


FORTRAN PROGRAMS FOR SIMULATING SOIL-MOISTURE 
AT THE SURFACE LAYER 
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1. PROGRAM TAYLOR. FORTRAN 
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OKieSfJAL PAeS JS 

OF POOR QUALITY 


C THIS PROGRAM GENERATES RAINSTORM EVENTS, STORM DURATIONS 
C AND INTERSTORM PERIODS WHICH PRESERVE THE HISTORICAL STATISTICS. 

C IT CALCULATES THE SOIL MOISTURE OVER A DEPTH CLOSE TO THE SURFACE 
C EVERY HALF HOUR .IT PLOTS THE EVAPOTRANSPIRATION 
C FUNCTION, THE SURFACE RUNOFF AND PERCOLATION FUNCTIONS. IT ALSO 
C PLOTS THE DAILY SOIL MOISTURE DURING THE RAINY SEASON LENGTH. 

C IT CALCULATES THE TOTAL STORAGE CHANGE, THE CUMULATIVE 
C EVAPORATION AND YIELD AT THE END OF EVERY RAINY OR 
C INTERSTORM PERIOD 

C IT HAS THE OPTION OF USING MANABE'S MODEL 
C TO CALCULATE THE MOISTURE FLUXES 
C THE VALUE OF THE POTENTIAL EVAPOTRANSPIRATION RATE 
C IS SET EQUAL TO ITS ANNUAL AVERAGE VALUE 

C CLIMATIC AND SOIL VARIABLES 

c epr«average annual evapotransplratlon rate(cm/day) 
c mtbwmean time between storms(days) 
c mtr^mean storm duratlon(days) 
c mpa»mean annual preclpitatlon(cm) 
c mtau*mesjn rainy season length(days) 
c ta«average annual air temperaturs(C) 
c mnu»mean number of storms per year 
c n^soll porosity 

c k1«saturated Intrinsic permeabll 1ty(cm2) 
c copore dlsconectedness Index 
c 2r»surfaco layer tk lckness(cm) 
c Mo“vegetat Ion cover 
c Kv»plant coefficient 

c koparameter of gamma dlstlbutod storm depths 
c Lamdaoparameter of gamma distributed storm depths 

real m1n,mo,m,n,nu,k1 ,mtb,mtr ,mh, In 

real sjk( 20) ,y 1 ( 20) , soj (20) ,a77(20) ,b77(20) ,b78(20) 

real da(365) ,SKP(&^"^) ,st(365) ,b79(20) ,a79(20) ,ys(20) ,yg(20) 

real a*^8(20) ,day(36£>^ 

f 1 1(d, ;o)»1./(d*(1.-'Sc)^*(1.45-.0375*d)+5./3. ) 
external plot_^$setup (descr Iptors) 
external plot^Sscale (descriptors) 
external plot* (descr Iptors) 
character * 10*xax 1 s , yax 1 s 
f 1(em)o10,**( .66+.55/em+. 14/em*^2.) 
kl 1»1 
ran«1 . 

print, 'To use Manabes parameterization type 2 , otherwise 1' 

Input, mnb 

If(mnb.eq.l) go to 3020 

print, 'Input the Initial soil moisture so' 

Input , so 
go to 3021 

3020 print, 'Input the average annual soil moisture so' 

Input, so 

pr1nt,'Input Time step (In days) ' 

Input , t Is 

C NR^Number of rainstorm events you want to generate 
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3021 print, 'Input NR' 

Input, NR 

print, 'Input storm properties k and Lamda' 

Input, xk, ami 
lf(mnb.Gq.2) go to 3040 

print, 'For dally fluxes ♦/pe 1,for half hour fluxes typo 2 ' 

Input, f 1 

print, 'To plot S(t) type 2|Otherwlse 1' 
input , lot 

print, 'For cumulative fluxes after each storm and interstorm period typo 2, otherwise 1' 
input ,ucu 

print, 'To plot S(t) for different values of 2r typo 2 .otherwise 1' 
input ,szr 

print, 'To print the cumulative fluxes only at the end of the rainy season ty, 5 e 2 , otherwise 1' 
input , feu 

3040 prInt.'To print thr* rainstorm events type 2 , otherwise 1' 
input , rae 
11*1 

3003 print,' epr , mtb , mtr , mpa , mtau , ta , mnu , n ' 
input , epr , mtb, mtr , mpa, mtau, ta.mnu.n 
if(mnb.oq.2) go to 3022 

2020 print, 'Mo, Kv.kI.c.Zr' 
i nput , VQ , vk , k 1 , C5 , zr 
if ( vg.oq. 1 ) stop 
1f(ran.oq,2) go to 3004 
if(dif.eq.2) go to 3004 

C J(s)»evapotranspi rat irn efficiency function 
C Ys(s)°surface runoff function 
C Yg(s)«ground water percollatlon function 

1000 print, 'To plot d(s) and y(s) typo 2 , otherwise 1' 

Input ,pl 

If(pl.eq.l) go to 3004 
if(kl 1.oq.2) go to 3004 

prInt.'To draw different curves for J(s) for different climates type 2, otherwise 1' 
input ,dl f 

double precision sum 1 .mean 1 ,mean2, means, 628 
double precision sum2 
double precision sum3 

3004 if (ran. eq. 2) go to 807 
3022 if(raG.eq.l) go to 42 

prlnt.'STORM DEPTH STORM DURATION TIME BETWEEN ' 
print,' (cm) (days) (days) ' 

42 il»1 

C ♦*♦♦•* *♦♦♦•***♦♦♦♦**♦**♦*•*•♦**♦♦*♦»♦**♦*♦**♦*♦♦*♦*<.♦♦♦♦*♦**♦♦ 

C GENERATION OF RAISTORM EVENTS 

C R1(I)astorm depth(cm) 

C R2(l)astorm durat lon(days) 

C R3( I )» interstorm durat ion(days) 

real R(500) ,WK( 1000) ,R1 (500) ,R2( 500) ,R3( 500) 
double precision DSEED 
DSEED-12376S.0D0 
A»xk 

Bo 1 . /ami 

call ggamr(DSEED,A,NR,WK.R) 
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do 5 l«1,NR 
R(I)«B*R(I) 

5 continue 
do 41 lal.NR 
R1(I)«R(I) 

41 continue 

DSEED-3478768.0D0 

A«1. 

B«mtr 

can ogamr(DSEED.A.NR,WK,R) 
do 7 I»1,NR 
R(I)°B*R(I) 

7 continue 
do 21 I°1,NR 
R2(I)t»R(l) 

21 continue 
DSEE0«649853.0D0 
A« 1 , 

B»mtb 

can ogamr(DSEED,A.NR.WK,R) 
do 9 I«1.NR 
R(I)°B-‘R(I) 

9 continue 
do 30 I»1 ,NR 
R3(I)-R(I) 

30 continue 
1f(ran.eq.2) go to 807 
If(rae.Gq.l) go to 3023 
go to 3024 

3023 1f(mnb.eq.2) go to 3025 
go to 807 

3024 do 11 I°1.NR 
write(6,17) R1(I),R2(I),R3(I) 

17 format(f 10.6. 4x,f 10.6, 4x.f 10.6) 
11 continue 

807 m'=2./(cs-3. ) 
d»cs“ 1 . /m- 1 
dE«2 . + 1 . /m 
f ledof ie{dE) 

c COMPUTE WATER CONSTANTS 


can WATCN(ta.sut .nu.gamsw) 

c COMPUTE CLIMATIC PARAMETERS 

del ta® 1 . /mtr 
mh»mpa/(mtou/(mtb+mtr) ) 
omnu*mtau/ ( mtb+mt r ) 
mi «mh/mtr 
eta* 1 . /mh 
olpna«'1 ./ml 
p1-3, 14159 
beta* 1 ./mtb 

c ♦♦*♦♦*•**♦♦*♦♦♦♦♦•*♦*♦*•«***•*♦*♦♦♦********♦«».*•♦♦*•••*•*••♦**♦* 

c COMPUTE DERIVATIVE OF J WITH RESPECT TO so 
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don-( 1.42B-0.037S*d) 

If (pi .oq. 1 ) go to 802 

koQ 

5O°0. 

809 4O«so+ 0.05 
go to 802 

802 ds"( 1 . -5o)**den 
ddG»da^d 
deno»dds+(5./3. ) 
denom»deno» ♦ ( 4 . /3 . ) 

5oo» 1 . -ao 

S0l*500**( -4 ,/3. ) 

dono3»2*8oo*deno 

dt°{2, 425-0. 0375<'d) 

so2»50o**dt 

dooa1»ao2*d*den 

nom« -denos -deos 1 

nom1»nom*so1 

dor onom 1 / ( donom* 3 ) 

f 1c«f Urn) 

si 1»sqrt(n/(k1*f 1c) )*sut/gamsw 
si 11"3l 1*8o**(-1 ./m) 
bk 1»k1 •gamsw/nu 

Slgc»n*otn*o2 . *bk 1 ♦si 1/(p1 ♦m*dal ta) ♦72000. 
6 1 gel -si gc**0. 3333333 
uars ig«H i gc 1 »dor 
s1a«5^n*bk1 • 88400* o1 l/(3*m*p1 ) 
s1gma"( s 1gc/dono»( 1.-so)^*2.)^*. 333333 
0»n1pho*bk 1 ♦86400* . 5^( 1 . +so^*cs) 
gl »a1ool0{8igma) 
xp«( 1 .766*g1 )+(0.980*(gl**2. )) 
xp1»- .806-xp 
CSI° 10. **xp1 
xp2«( 1 .96*gl )+l .766 
U»-dGrs Ig*xp2/s igma 
co*ft1pha*86400*bk1/2. ♦cs*so**(cs- 1 . ) 
co1»U-co 

C2>*co1 ♦CSI ♦oxpC -g) 

C38“mtou*86400*bk 1 ♦cs/mpa*so**(ca- 1 . ) 
C3»C.;.8/2. 

If(vg.oq.O) go to 80 
go to 90 


GO E*2 . ♦bota*n*bk 1 *3 1 1 ♦f 1ed/(p1 *m*epr*^2. )*86400*80**(d+2. ) 
If (E. go. 88, ) E«88. ' 

z 1«( 1 .+E*sqrt(2 . ) )*exp( -E) 
z2«gomma( 1 ,9)-gamt( 1 ,5,E) 
z2-z2*8qrt(2. *E) 
sj n 1 . -z1+z2 
If (pi .eq, 1) go to 803 
k»k+1 
sjk(k)«sj 

1f(k.Qq,20) go to 804 
go to 805 

803 ag»gammo( 1 .5)-gamt( 1 ,5,E) 
g1»oxp( -E)*sqrt(2. ) 
g2“E*8qpt(2. 
g 2 «g 2 * 0 xp( -E ) 
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quality 


g3"O0»sqrt(2, )/(2. *»qrt(E ) ) 

04»oxp( -E)«5qrt(E)*i5qrt(2*E) 

0g»-QUQ2+g3-g4 

E 1 1»2* *beta*n*bkl*9l 1*f l0d/(p1*m*apro»2. ) *80400 
E12»(d+2.)*oo**{d4.i.) 
d«r1 J»gg*E11*E12 
C1“dor < J 

if(Cl.lo.O) C1»0.0 
go to 100 

00 B«( . -vg)/( 1 . + (vg*vk) ) 

B“B+(vk*vg**2. )/(2.*(1.+(vg*vk))**2. ) 
C»1*/(2.*(vg*vk)**2.) 

E 1«2 . *bota*n*bk1 ♦si 1 *f iod/(p1 *m*epr**2. ) *86400 

E»2 . *bQta*n*bk1*5l 1*f ied/(p1 *m*epr**2, )*86400*aa**(d*»-2. ) 

o1oB*((vg*vk)+1) 

o1»-oUsqrt(B*2. ) 

o1 1«B*E*sqrt(2. *B) 

o1»o1-o1 1 

o1 »o1 *QXp( -B*E ) 

o1“ot*E1*(d+2. ) 

o1»o1*(so**(d+l . )) 

o2o-vg*vk*C 

o2»o2+sqrt(2*C) 

o2oo2-(C*sqrt(2*C)*E ) 

CSS«G*£ 

1f(C88.ge.88) C88»88. 
o2°o2*Gxp(*C88)*E1*(d'*-2. ) 

02»o2*(oo#*(d+1 . ) ) 

CE°0*E 

BE»D*E 

o1«(vg*vk)+ 1 . 
a2»E*Qqrt{2. *B) 
o3ao1+a2 

if (BE.ge.88. ) BE«88. 
a3<>a3«oxp( -BE ) 
a4»V0*vk 

o4»a4+(E*sqrt(2. *C) ) 

If (CE.ge.88. ) CE«88. 
a4»o4*exp( -CE) 

a5»gamt( 1 . 5 ,CE ) -gamt ( 1 . 5 .BE ) 
o5«a5*sqrt(2. *E ) 
a6«a3-G4 -aQ 

a6»a6*(1. 'Vg)/( 1 . -vg+(vij>vk) ) 
sj* 1 . -aG 

If (pi .oq. 1 ) go to 806 

k»k+1 

sJU(k)»sj 

if(k.eq.20) go to 804 
go to 805 

806 o3cgamt( 1 .5,CE)-gamt( 1 ,5,BE) 
o3»o3*sqrt(2.*E1) 

03-«( 1 , + d/2. )*o3*( so**(d/2 . )) 
o31--C*E1*(so**(d^2. )) 

031=^(C<i* 1 .5)*©xp(o31 ) 
o32«-B*E1*(so**(d+2. )) 
o32“(B**1 .5)*exp(o32) 
o33»o31-o32 
033»033*(E1**1,5) 
o33«o33*(2.+d) 
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o33»o33*(bo**( ( 1 .6*c0'»-2. ) ) 

o33«o33*oqrt{2. ♦£) 

o3*o3+o33 

derj oQl -o2-o3 

dor,1 BdorJ ♦ ( 1 . -vg) 

darj »-dopJ/( (vQ*vk)+1 , -vg) 

Clodorj 

if(Cl.lo.O) C1-0.0 
B28»mtau*bk i •86400/mpa^so**C5 

r :ioDerWative of 0 with respect to s 
c C2«0er ivative of Yo with respect to s 
C C38°Derivativ© of Yg with respect to s 
C 0joiJ(5O) 

C si11“psi evaluated at so 

C bk1«3oturatod hydraulic conductivity (cm/aec) 


100 print 101 .Cl ,C2.C38,6j ,5i ll.bkl 

101 format(3hC1-,f 10.6.4x.3hC2»,f 10.6,4x,3hC3«,f 10.6,4x,2hd-,f 10,6.4x,3hMH«‘.f 10.2,4x,f20. 10) 
SK»so 

B 04 p»fflpa/(mnu*mtr ) 

C1«C1 -opr 
B1«sj -epr 

if (pi .oq. 1 ) go to 808 
80 * 0 . 
k«0 

811 so«so+0.05 
ds*d*( 1 -so)»-d0n 
dono*ds+(5,/3. ) 

sigma»(s1gc/dono*( 1 . -so)-^2. )-*. 333333 

808 B22*sigma*-(-s1gma) 
sigmoslgma+1 . 

B22«B22*gamma(sigm) 

B2*B22*ex;j( -g-(2*s1gma) ) 

B28 *m t au *bk1 -86400/ mpa* 50* -cs 
B4"B2*p 

B5*B28*p*ninu-mtr/mtau 

If(pl.eq.l) go to 809 

k»k+ 1 

ys(k)°B4 

yg(k)*B5 

soj (k)*so 

1f(k.eq.20) go to 810 
go to 811 

809 1f(ucu.oq.2) go to 1816 

print.' S(t) 1(cm/day) Et(cm/day) y 1old(cm/day ) DAY ' 

go to 1815 

C 414 *«««**#* 

C CALCULATE THE SOIL MOISTURE CONCENTRATION AND THE 
C CUMULATIVE EVAPORATION AND YIELD AT THE END OF 
C EVERY RAISTORM AND INTERSTORM PERIOD 

C ****«**«^#*««#«*»4.4.*«i«*«****«**-»«*4i**«*«*4<**«**«i»«>****il>«****« 


1816 print, 'SOIL, MOIST. CUM.EVAP. CUM. YIELD 
1815 If (pi .oq. 1 ) go to 812 
810 1f(kl1.oq.2) go to 3001 
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Of* POOH 


c •♦♦*♦♦♦♦•♦*♦♦♦♦ 

C Plot vl VinVaUS tt 


c«n plot $8«tup(' *• 'ti' , 'US 1 .0.0,0) 
can plot $»e«il«iiO. , 1 . ,0. . 1 . ) 

3001 IK) 

HI 1«3 

do « 1 a J « 1 , 20 

b7?(0**»JU(J) 
ttr' M 1 )*»»0,1 ( J ) 

D13 con t1 rum 

eoM plot^ (fl7M>7?,ao. 1 . ' M 
1f(cHf.«q/l) QO to 3002 
00 to 3003 

3002 rond(».) 


C PLOT Yo AND Vg VtPSUS « 
C 


con plot $»otup(' S'SOU. MOISTURES 'SUPPACE RUNOPP M ,0,0.0) 

Cftll plot $ftonlo(0. , 1 . ,0. ,2. ) 

<«0 

do 6 In j « 1 , 20 

1 - H 1 

fc)7(H \ )«V8( J ) 

«7B( 1 )»*(Soj( J ) 

614 con 1 1mm 

Cftll plot.(«70.t)70,20. 1,' •) 

rmnd(n. ) 

cnll plot^1.8Qtup( * S*S0U M01STURI!S ^aROUNDWATTR RUNOHP M ,0,0,0) 

cni 1 plot $ncnU(0. , S ,0. ,2 ) 

1«0 

do 034 j«i,ao 

1 « 1 M 

b73U)«yg( J) 
ft73( 1 )>*ttOj ( J ) 

634 cent 1 mm 

G«l I plot («7R,b70,ao. 1, ^ ' ) 
go to 1000 

612 ir'last' oq. 1 ) go to 617 
do 2001 11 **1.2 
prMnt , ' Input /iMcin) * 

Input ,xr' 

6 <7 o«n*jer' 

Ot«t 18 

K-0 

KP*0 

I«0 

IM«0 


SK3-0.0 
SK2**0.0 
LMM-0 
y 1oldc«0.0 
ovnpc*»0. 0 

400 1 f ( ucu . oq . 1 ) go to 401 
1 f ( 8 / 1 ' . oq . 2 ) go to 401 
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if(fcu.eq.2) go to 401 
wr 1 te(6, 1701 )SK,ovapc,y ieldc 
1701 format(f8.5,4x,f8.5,4x,f8.5) 

C CALCULATE THE VALUE OF SOIL MOISTURE EVERY HALF HOUR 
C DURING A PRECIPITATION EVENT 


401 Dt1»0. 
yt»0.0 

sia1°sia*f 1 i (d,SK) 

B 1a2«2*( 1 . -SK)*sqrt(s1a1 ) 

Ao»bk1 ♦86400/2. 

If(SK.le.O) go to 216 
ao1»Ao*( 1 .•♦•(SK**cs) ) 
go to 216 

215 ao1*Ao 

216 I«I+1 
r2-R2(I) 
in«R1(I)/r2 

Tol »2* 1n*( In-aol ) 
to2°s1a2+*2./To1 
to3«2. ♦( In-aol ) 
to4^ 1 .+{eo1/to3) 

Toato2*to4 

300 Dt1«Dt1+Dt 

If (Dtl ,ge.r2) go to 200 

LM-LM+I 

1 f (Dt 1 .ge.To) yt“1 
if (SK2. 1 t.SK3) yt«0.0 

SKI ^SK+C 1n-p*( (B2*yt ) + (B28*mnu*mtr/mtau) )-p*(SK-so)*( (C2*yt ) + (C3*mnu*mtr/mtau) ))*Dt/a 

SK2*5K1 

SK3=SK 

If (SKI .ge. 0.999) go to 211 
go to 212 

211 SK1»0.999 
y lei d* in 

y ieldc«y1e1dc+( In^t Is) 
go to 213 

212 y1eld-p*( (B2*yt) + (B28=<‘ninu*mtr/fntau))+p<'(SK-so)*((C2*yt) + (C3*mnu*mtr/(ntau)) 
y lei dc«y iel dc+(y iel d>*»t 1 s) 

213 SK»SK1 

if(fl.eq.l) go to 250 
1f(szr.eq.2) go to 300 
write(6,210) SK, in, yield 
210 format(f8.5.4x,f8.5,22x,f8.5) 
go to 300 

250 t1ss«1./tis 

i f (LM.ge. t iss) go to 251 
go to 300 

251 LM»0 
LMM«LMM+1 

if (LMM.gt .mtau) go to 900 

KP»KP+1 

SKP(KP)=*SK 

da(KP)»LMM 

1f(ucu.eq.2) go to 300 
1 f (szr ,eq. 2) go to 300 
write(6,252) SK, in, yield. LMM 
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252 format(fB.5,4x,f8.5,22x,fB.5.9x, 15) 
go to 300 

200 1f(ucu*eq,1) go to 201 
1f(fcu.eq.2) go to 201 

wri to(6, 1700) SK.yleldc.yt 
1700 format(f8.5,‘16J<,f8.5,4x,f3. 1) 

C CALCULATE THE VALUE OF SOIL MOISTURE EVERY HALF HOUR 
C DURING AN INTERSTORM PERIOD 

201 Dt1«0. 

500 Dt1»0t1+Dt 
r3°R3(I ) 

If (Dt1 ,ge.r3) go to 400 
LM»LM+1 

evap=»B1 + (C1*(SK-so)) 

if (evap.ge,epr) go to €00 

evapps»evap/epr 

If (evapp. le. vg) go to 701 

SK1»SK-(evap4‘(B28*p*mnu*mtr/mtau)+(C3*p*mnu*mtr*(SK-so)/mtau) )*Dt/a 

evapc«evapc+ ( evap* 1 1 s ) 

go to 700 

600 evap»epr 

evapc=evapc+ ( evap • t 1 s ) 

SKl»SK-(epr*Dt/a) -( (B28^p^mnu*mtr/mtau)+(C3*p*mnu*mtr*(SK-so)/mtau) )*Dt/a 

go to 700 

701 evap»8pr^vg 

evapc=evapc+(evap*t 1 s) 

SK 1°SK-(evap*Dt/a)-( (B28*p*mnu*mtr/mtau)+(C3*p*mnu^mtr»(SK"So)/mtau) )*Dt/a 
700 y 1eld°(B28*p»mnu*mtr/mtau)+(C3*p*mnu*mtr*(SK“So)/mtau) 

If (yield. le. 0.0000001) y lei d»0. 0000001 
y 1eldc*y 1eldc+(y leld^t Is) 

SK«SK1 

If (f 1 .oq. 1 ) go to 750 
1f(szr.eq.2) go to 757 
wr 1 te(6 , 220) SK , evap , y lei d 
220 format(f8.5, 16x,f 8 .5, lOx, f8.5) 

757 K«K+1 

if (K.ge. 1000) stop 
go to 500 

750 t1ss=1./t1s 

If (LM.ge.tiss) go to 751 
go to 500 

751 LM=0 
LMM»LMM+1 

If (LMM. le.mtau) go to 901 

write (6, 905) SK.evapc.y leldc 

905 format(f8.5.4x.f8.5.4x,f8.5) 

go to 900 

901 KP«KP+1 

SKP(KP)=SK 

da(KP)«LMM 

1f(ucu.eq.2) go to 500 
1f(szr.eq.2) go to 500 
wr1te(6,752) SK, evap, yield, LMM 

752 format(f8.5, lex.fS.S. 10x,f8.5,9x, 15) 
go to 500 

900 1f(szr.eq.2) go to 2008 
1f(ran.eq.2) go to 2031 
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C CAUCJLATE THE STATISTICAL PROPERTIES OF THE GENERATED 
C RAINSTORM EVENTS 

C <»****#i****4»********»#*«#M4n«i**i*i*»*i«i*«i*r**«4i*4t*>#**^*****Hr*******m# 

print , 'Statistical properties of the simulated rainstorm characterlst ics' 
sum1«O.ODO 
sum2*0.0D0 
sum3»O.ODO 
do 1001 IL«1,I 
suml«sum1+Ri(lL) 
sum2«sum2'fR2( IL) 
sum3»sum3+R3( IL) 

1001 continue 

mean1 "suml/Cf loat ( I ) ) 

mean2«sum2/(f loat( I ) ) 

mean3>sum3/ (float(I)) 

var 1 »0.0 

var2«0.0 

var3=0.0 

do 1002 IL°1.I 

var 1 =var 1 + ( (R1 ( lL)-mean1 )>»"»'2. ) 
var2“var2+((R2(lL)-mean2)**2. ) 
var3»var3+( (R3(IL) -means )*+2. ) 

1002 continue 

var i 1=var1/f loat{ I- 1 ) 
var i 2“ var 2/ float ( I - 1 ) 
var 13®var3/f loat( I - 1 ) 

print, 'AVER. h(cm) AVER . tr (days) AVER . tb(days) ' 

wri te(6, 1003) mean 1 , mean2 , mcan3 

1003 format(f 10.6, 6x,f 10.6, 6x.f 10.6) 

print,' VAR.h VAR.tr VAR.tb' 

print 1004 , var i 1 .var 12 , var i 3 

1004 format(f8.2,4x,f8.2, 10x.f8.2) 
ran*2. 

2031 if(lot.eq.2) go to 2030 
go to 2020 
2030 read(5.) 

2008 if(il.gt.l) go to 2003 

C PLOT THE SOIL MOISTURE CONCENTRATION WITHIN THE 
C LAYER OF THICKNESS Zr VERSUS TIME DURING THE 
C RAINY SEASON LENGTH 

call plot_$setup( ' ', 'DAYS' , 'SOIL MOISTURE 1 ,0,0,0) 
call pi ot_$scal e( 1 , , 220. ,0, , 1 . ) 

2003 ino 

do 910 j»1,LMM 
i = i + 1 

st(1)=SKP(J) 
day( i )»daO ) 

910 continue 

if ( il .eq. 1) go to 2004 

if(il.eq.2) go to 2005 

2005 call plot_(day , St ,mtau, 3, ' . ' ) 

go to 2001 

2004 call plot_(day, st.mtau, 1, ' ') 
if(szr.eq.l) go to 2000 
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2001 continue 

C CALCULATE THE MOISTURE FLUXES USING MANABE'S PARAMETERIZATION 

Q ^Hlt^^(H¥^^*^*^4‘***^**^^*^l^***^***^***^***‘^*******^*********>***'^*^<¥******* 

3025 If(mnb.eq.l) go to 2000 

print, 'S(t) CUM.EVAP. CUM. YIELD' 

SK“SO 
Dt«1 ,/48. 

I«0 

y 1eldc*0. 0 
evapc«0.0 
Dtl 1«0.0 

3031 wr1tG(6,3033) SK , ovapc, y lei dc 
3033 format(f8.5.4x.f8.5,4x,f8.5) 

Dt1«0.0 

1 = 1 + 1 
r2=R2(I) 

1n=R1(I)/r2 

3028 Dt1»Dt1+Dt 
Dt 1 1»Dt 1 1+Dt 

If (SK.ge.0.42) go to 3029 
SK1=SK+1n^0t/(n*100) 

SK«SK1 

If (Dtl .ge.r2) go to 3027 
go to 3028 

3029 y ield=( 1n-epr)*Dt 
y 1 e 1 dc=y 1 el dc+y 1 e 1 d 

If (Dtl .ge.r2) go to 3027 
go to 3028 

3027 wr 1 te(6 , 3030) SK, Gvapc,y lei dc 

3030 forfnat(f8.5.4x,f8.t;,4x.f8.5) 

Dt1=0.0 

r3=R3(I) 

3032 Dt1«*0t1+0t 
Dtl 1»Dt11+Dt 
evap»epr 

if (SK. It. 0.3 15) evap»epr>*SK/0.315 
SKI^SK-evap^Dt/Cn* 100) 
evapc=evapc+(evap*Dt ) 

SK=SK1 

If (Ot 1 1 .ge.mtau) stop 

If (Dt1.ge.r3) go to 3031 
go to 3032 
2000 read (5.) 

stop 

end 


subrout 1 ne WATCN( ta , su . nu , gamsw ) 


C ♦****♦♦♦****♦♦♦♦-.•♦*♦*♦♦***♦♦*♦♦*♦♦**♦*♦*♦♦***♦♦•♦ *♦*l♦♦*♦**♦♦♦* 

real nu.nut 

dimension sutt( 1 1 ) ,nut( 1 1 ) ,gamst( 1 1 ) 

data sutt/75.6,74.9,74.2.73.5.72.80,72. 1.71.4,70.7.70.0,69.3.68.6/ 
data nut/17. 93e-3, 15. 18e-3. 13.09e-3, 11.44e-3, 10.086-3,8.940-3, 

& 8. e-3, 7. 2e-3. 6. 530-3,5. 976-3,5.946-3/ 

data gamst/O. 99987, 0.99999999, 0.99973, 0.999 13, 0.99823, 0.99708, 
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& 0. 99568,0. 99406, 0.9922&. 0.99025, 0.98807/ 

If ( ta.gt.50. )go to 10 
1 ta» If 1x( ta* . 2)+1 
f rac»ta- float (5*( 1 ta- 1 ) ) 

1 tal •• 1 ta+ 1 

8ut«(sutt( 1 ta 1 )“5utt( 1 ta) )*0. 2*froc+sutt( 1 ta) 
nu«(nut( 1 tal )-nut( 1 to) )*0. 2*f rac+nut( 1 ta) 
gamsw"( (gamstC 1 tal )-gamst( 1 to) )* . 2«f rac+gomstC 1 ta) )*980, 
return 

10 sut«autt( 1 1 ) 
nu-nut( 1 1 ) 
gomsw«gamst( 1 1 ) 
return 
end 

c this function computes the gamma Incomplete function 

function gamt(a.x) 

1f(x.oq.0)go to 40 
If (x.gt . 100)go to 50 
Bum« 1 . /a 
ana 1 .0 
old«sum 

33 ol d«old*x/(a+on) 

If (old/sum- 1 . 0 - 6 ) 20 , 10,10 
10 on«an+ 1 . 
sum-sum+ol d 
If (an-300. )33.33. 12 
12 continue 

20 xxx*(a*alog(x)+alog(sum)-x) 

1 f ( XXX . 1 t ■ -80. )go to 40 

gamt"(oxp(xxx)) 

go to 60 

40 gamt»0.0 

go to 60 

50 Qamt»gamma(a) 

60 return 
end 

c This function computes the gamma function by a Stirling approx. 

function gamma(y) 
x«y+ 1 . 
p1«3. 14159 
st1r1«1./(12.*x) 

St1r2°1./(288.*x**2.) 
st 1r3«-139./(51R-*0.*x**3. ) 

6t1r4--571 ./( 2488320. ♦x**4. ) 

St 1 r«« 1+st 1 r 1+st 1r2+st 1 r3+st 1r4 

gamma«exp( -x)*x**(x- .5)*sqrt (2 . ♦pi )*st Ir/y 

end 

function f1e(d) 
dimension y(6) 

data y/0. 18,0. 11,0.077,0.056,0.044,0.034/ 

1 f (d. gt . 7 . )go to 10 

X"d- 1 . 

1»1f 1x(x) 
f rac«x-f loat( 1 ) 
y1-alog(y( 1 )) 
y2«alog(y ( 1+ 1 ) ) 
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f iQ«exp( (y2-y1 )*f rac+y 1 ) 

return 

10 f1e».034 

return 

end 
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13 

r ♦’.* V * 
^ 0 ^ 


C THXS PROGRAM CALCULATES THE AVERAGE SOIL MOISTURE CONCENTRATION 
C OVER 1m DEPTH EVERY HALF HOUR , DURING AN EVAPORATION PERIOD WITH 
C A CHANGING VALUE OF THE POTENTIAL EVAPORATION RATE. 

C IT ALSO CALCULATES THE SURFACE TEMPERATURE USING THE FORCE- 
C RESTORE METHOD OR THE THERMODYNAMIC EQUILIBRIUM EQUATION. 

C THE POTENTIAL EVAPORATION RATE IS CALCULATED EITHER USING 
C PENMAN'S EQUATION OR THE AERODYNAMIC EQUATION. 

C ATMOSPHERIC INSTABILITY CRITERIA ARE USED. 

C THE SHORT-TERM WATER AND THERMAL BALANCES CAN BE SOLVED SIMULTANEOUSLY 
C AND THE SOIL MOISTURE CONCENTRATION AND SURFACE TEMPERATURE 
C CAN BE CALCULATED AND PLOTED 

C THIS PROGRAM READS FROM FILE 31 THE METEOROLOGIC 
C VARIABLES AND THE SOIL MOISTURE CONCENTRATION MEASUREMENTS, 

C WHICH ARE GIVEN EVERY HALF HOUR. 


C THE CLIMATIC VARIABLES AND SOIL PARAMETERS USED AS INPUTS 
C TO THIS MODEL ARE DESCRIBED BELOW. 

c epr«annual average potential evaporation rato(cm/day ) 
c mpa*«mean annual preclpitation(cm) 
c mtr=mean storm durat lon(days) 
c mtau^iflean rainy season length(days) 
c mnu»mean number of storms per year 
c Jaevapotranspirat Ion efficiency 
c Cloderivative of J with respect to s 
c C3»der Ivatl ve of percolation rate with respect to s 
c so«average annual soil moisture 
c SK^Inltlal soil moisture at 1m depth 
c naporosity 

c Zr«surface layer thichness 

c K( 1 )osaturated hydraulic conductivity (cm/sec) 
c capore d1 sconectedness Index 
c a(1,1)“net radlat l 9 n( 1 y/mln) 
c a(1.2)»a1r temperature(C) 
c a(1,3)awater vapor pressure of air (mb) 
c a(1,4)"w1nd speed(cm/sec) 

c a( 1 ,5)«average soil moisture content In 0-10cm 
c a( 1 ,6)«average soil moisture content In 10-50cm 
c a( 1 , 7)aaverago soil moisture content In 50- 100cm 
c a0.8)»ground temperature at 1cm (C) 

c Tg»ca1culated surface temperaturo(C) 
c T2“deep soil temperature(C) 

c (cH)n»drag coefficient under neutral conditions 

real a ( 337 , B ) . ep 1 ( 337 ) , hr ( 337 ) , ASK ( 337 ) 
real epp(337) ,hrr(337 ) 

real AsK( 337 ) , SK2( 337 ) . hr 1 ( 337 ) , SK3 ( 337 ) . hr2( 337 ) 
rea 1 AsK 1(337), TgCC ( 337 ) , TgC 1 ( 337 ) , TgKK ( 337 ) 
real TgCM(337) 

external plot_$setup (descr Iptors) 
external plot_$scale (descriptors) 
external plot~( descriptors ) 
real mpa,mtau,mtr ,mnu,k1 ,m 
f1(em)-10.**( .66+.55/em+, 14/em»*2. ) 

double precision B ,ga, gd,gd1 ,bet ,deno,T,es,d1f ,H,nom,ep,B28 
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print, 'epr ,mpa,mtr, fntau,mnu,ii, Cl ,C3,so,SK,n,2r,K( 1 ) ,c' 

Input ,epr,mpa,mtr,mtau,mnu,sj ,C1 ,C3 .ao.SK,un, 2 P,bk 1 ,cs 

roadOl, ) ( (a( i . j ) , i * 1 .337) . J « 1 ,8) 

print, 'To print f11e31 typo2, otherwise type 1' 

Input , ty 

print, 'To plot ep typo 1. otherwise type 3' 
input ,pr 

prInt.'To print ep write 2 .otherwise 1' 

Input ,pr 1 

prInt.'To print soil moisture type 2 .otherwisel' 

Input ,pt 

print. To calculate the surface temperature type 2, otherwise 1' 

Input , tmr 

If(tmr.eq.l) go to 200 

print.'input the Initial surface temperature Tg, T2 (In degrees Colclua) and (cH)n' 
Input, TgC,T2,cHn 

prInt.'To solve simultaneously the equations for soil moisture 
& and temperature using the aerodynamic equation and 
& the Instability criteria type 2 . otherwise 1' 

Input ,aer 

prInt.'To use the thermodynamic equation type 2 , otherwise 1' 
input , thm 

if(thm.Qq.l) go to 305 
print.'input k(1),Ta' 
input ,k1 , ta 

m*2 ./(cs-3 . ) 
f 1c»f 1 (m) 

C COMPUTE WATER CONSTANTS 

call WATCN(ta,sut,nu,gamsw) 
si 1»sqrt(un/(k1*f 1c) )*sut/gamsw 
305 T2»T2+273.16 
TgK°TgC+273. 16 
TgFB(9.*TgC/5. )+32. 

TgKK(1)»TgK 

200 If(ty.eq.l) go to 41 
do 40 1-1,337 

wrlte(6,20) a(1,1),a(1,2),a(l,3),a(l,4),a(1,5),a(1,6).a(1,7),a(1,8) 

20 format(f 10.4,2x.f 10.4.2x,f 10.4,2x.f 10.4.2x,f 10.4,2x,f 10.4,2x,f 10.4.2x,f 10.4) 

40 continue 

41 h»-0.5 
sum-0.0 

do 46 1-1,337 

C CALCULATE ep USING PENMAN"S EQUATION 

ga«(a( 1 ,2)*0.013)+0.42 

gd 1 - 1 , /ga 

gd«gd1-1. 

dano“597. *gd1 

bet»200./0.03 

bet«a1og(bet) 

bet=bet**2. 

B-10.**(-7) 

B«1 .222*B 
B-B*a( 1 ,4)*60. 
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B»B/bet 

T-273, ie+a( 1.2) 

68-273, 16/T 
es-1 .-08 

©5-08*(5. 006504-19. 83923) 

05-exp(es) 

e8l«(273. 16/T)**5, 00650 
es-its-esl-C. 1 1 
d1f'>ea-a( 1 ,3) 

H°597.*B-cMf#gd 
nom-H+a( 1,1) 

64 op-nom/deno 
ep-6p»60. *24. 
h-h+0. 5 
epK 1 )°ep 
hr( 1 )-h 
sum-sum+epi ( 1 ) 

1f(prl.0q.1) go to 46 
wr1te(6,45) ep 

45 formatlf 10.4) 

46 continue 
avop-8um/337. 
write (6. 60) avep 

60 format(2x, f 10.4) 

1f(pr.eq,2) go to 61 

cal 1 plot_$30tup( 'Potential Evaporation' . 'Hours' , 'ep' , 1 ,0,0.0) 
cal 1 plot3scale(0, . 168. ,-0.15,2.) 

61 1»0 

do 51 jni,337 
1-1 + 1 

epp(l)-epUj) 
hrr( 1 )»hr(. ) 

51 continue 
1f(pr.eq,2) go to 63 
call plot_ (hrr ,epp, 337 , 1 , ' ') 

C CALCULATE THE UPDATED SOIL MOISTURE 

63 a1-un*zr 
hr8»-0, 5 
p-mpa/(mnu*mtr ) 

B28»mtau*bk1*86400./mpa*so**cs 
Dt-1 ./48. 

C33-C3/2. 
do 100 1-1,337 

ASK(1)»(0. 40-a(1,5)) + (0.40*aCl,6)) + (0.50*a(1.7)) 

S1-3j-epr 

c1»C1*epr 

yt"0.0 

evap«B1+(c1*(SK-so) ) 

1f(1.1e.196) go to 108 
1f(evap. It.epl(D) go to 107 
108 evap-epl(l) 
yt«1 .0 

107 AsK(1 )-ASK( 1 )/0.35 
ASKK-^ sK(1) 

SK2(11»SK 

hr1(1)^Hr8+0.5 

SK1-SK-(evap+(B28*p*mnu*mtr/mtau)+(C33-p*mnu*mtr*(SK-so)/mtau) )*0t/a1 
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if (pt .Qq. 1 ) go to 102 
wr1to(6, 101) ASKK.SK.yt 

101 format(f 10.4,4x,f 10.4,4x,f3. 1) 

102 SK*-SK1 
hr8«»hPl( 1 ) 

100 continue 
1f(tmr.eq.2) go to 290 

C PLOT THE CALCULATED AND MEASURED SOIL MOISTURE CONCENTRATION 
C AT DEPTH OF U» . USING PENMAN"S EQUATION TO CALCULATE op 

C «i<i*****««*«*4i*«*««i«*«*****«w****i*«****i»**«»****»****»***4<«*«*«** 

con plot $setup(' '.'HOURS', 'SOIL MOISTURE' , 1 ,0,0,0) 
cal 1 plot“$sca1a(0. , 170. ,0.62,0.70) 

1-0 

do 110 J-1,337 
1-1 + 1 

AsK1(1)-AbK(J) 

SK3( 1 )-SK2(J ) 
hr2( 1 )-hr1(J ) 

110 continue 

call plot (hr2.SK3,337,3. ' . ' ) 
call plot" (hp2.AsK1.337.1.' ') 

290 0t»1800. 
if (tmr.eq. 1 ) go to 210 

Ev1“- 10. 

SUM-0.0 
L- 1 

print, 'Average Dally Evaporation Rate(cm/day) ' 

SK-SK2(1) 

C *«*****************************i**** 

C CALCULATE op USING THE AERODYNAMIC EQUATION AND THE 
C ATMOSPHERIC INSTABILITY CRITERIA( surface roughness 0.05cm) 

C 


do 2S0 1-1,337 
TgA-a( i ,2)+273. 16 
SK2( 1 )-SK 

est-6. 11+(0.339*(TgF-32. )) 
if(1.1e.196) go to 260 
Ev11-B1+(cl-(SK2(l)-so)) 

Evl-Evi 1/86400 

260 Ri-2.-981*100.*(TgA-TgK)/((TgA+TgK)»(a( 1 ,4)**2.)) 
If (R1 .go. 0.2) rat-0.0 

if (Ri . 1 t,0.2.and.Rl .ge.O. 1) rrt-(-2. *R1 )+0.4 
if (Ri . 1 1 .0. i .and.Rl .ge.0.0) rat-(-8.*R1)+1. 
if (Ri . It.O.O.and.Ri .ge.-O. 1) rat- 1.30 
1f(Ri. It. -O. 1. and.Rl .ge.-0.2) rat-1.8 
If (Ri . It. -0.2.and.Rl .ge.~0.3) rat-2.2 
If (Ri . lt.-0.3.and.Rl .go. -0.4) rat-2.45 
if(Rl .U.-0.4) rat-2. 7 
cH-cHn+rat 

Ev-cH-(730.5e-9)*a( 1 ,4)-(est-a( 1,3)) 

1f( 1 . 1e. 196) go to 262 
If(Ev.ge.Evl) Ev-Evi 
if(thffl.eq.l) go to 262 
If(Ev.lt.Evl) go to 262 
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CALCULATE THE SURFACE TEMPERATURE USING THE 
THERMODYNAMIC EQUILIBRIUM EQUATION 


O8t1«Ev1/(cH*(730.5o*9)*a( 1 ,4)) 

08t«e»t1+a( J ,3) 

8m»8i1*SK**(-1./m) 

OX-98 1 . *8 1 1 1 / ( ( 2 . 8760+6 ) *ToK ) 

rh-Qxp(ex) 

TgF1«e8t-(rh*6. 11)+(0. 339*32. *rh) 

TgF«TgFi/(0.339*rh) 

J-1+1 

TgKK(J)-(5.*(TgF-32.)/9. )+273.1G 
TgK-TgKK(J) 

262 if(aer.cq.l) go to 261 

SK1«SK-( (Ev*064OO. ) + (B28*p*mnu*nitr/n)tau)+(C33*p*mnu*mtr*(SK-ao)/mtau))/(48. *a1 ) 
SK^SK 1 

If(Ev.U.Evi) go to 261 

1f(1. 10.196) go to 261 

If (thtn.oq.2) go to 250 

261 HS-cH*(285.48o-6)*a(1 .4)*(TgK»TgA) 

Le“597.3-(0.57*TgC) 

G-(b( 1 , i )/60. )-Hs-( Le*Ev) 

J-1 + 1 

C COMPUTE SURFACT TEMPERATURE 

TgKK(J)-TgKKH ',^ (2. •G*0t/7.37)-((72.72e-6)*Dt*(TgKK( 1 )-T2)) 

TgC-TgKK(J 16 

TgF»(9.*TgC/5. )+32. 

ToK-TgKK(J) 

T2-T2+(G*Dt/249.68) 

SUM-SUM+(EV*8S400. ) 

L-L+1 

If (L. It. 48) go to 250 
AEv-SUM/48. 
wr1to(6.4C0) AEV 
400 format(2x,f8.4) 

L-1 

SUM-0.0 
250 continue 

c *..******.*••*•**•**•*.•.•••••*••*•••*••**•*•*••.•*•*•••••••••* 

C PLOT CALCULATED AND MEASURED SURFACE TEMPERATURE 
C USING THE AERODYNAMIC EQUATION 

call ptot_$setup( ' 'HOURS' . 'SURFACE TEMPERATURE 1 ,0,0,0) 
ca] 1 plot""$scale(0. , 170- , -2. ,40. ) 
do 270 1»iT337 
TgCC(i)«TQKK(1)-273.16 
270 continue 
1 »0 

do 280 J« 1,337 
1«1 + 1 

TgCM(1)=a(J,8) 

ToC1(1)»TgCC(J) 
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hp2( 1 )«hr1(J ) 

2BO continua 

call plot (hr2,ToC1 ,337,3, ) 
call plot_ (hr2,T0CM, 337,1,' ') 

(f(aar.oq.T) go to 210 

read(S, ) 

call plot $S8tup(' '.'HOURS', 'SOIL MOISTURE ', 1 ,0,0,0) 
cal 1 plot“$scalQ(0. . 170. .0.61,0.70) 

1-0 

do 300 J-1 ,337 

1-1 + 1 

SK3(1)-SK2(J ) 

AsKKi )-A5K( j ) 
hr2(1)-hr1(j) 

300 continue 

C PLOT CALCULATED AND MEASURED SOIL MOISTURE DERIVED BY 
C USING THE AERODYNAMIC EQUATION FOR ESTIMATING op 
C **4i«if***-**')i*-**-*«*****-***4i*««*«»***4<**4i**«****ili4i-*iti-«-«i****»* 

call plot (hP2,SK3,337,3,'.') 
call ploC (hr2,ABK1, 337,1 ,' ') 

210 stop ” 
end 

GUbhout Ine WATCN( ta, sut ,nu,samsw) 
real nu.nut 

dimension sutt( 1 1 ) ,nut( 1 1 ) ,oamst( 1 1 ) 

data autt/75.6,74.9,74,2,73.G.72.80,72. 1 .71.4,70.7,70.0,69.3,68.6/ 
data nut/ 17.930- 3, 15. 18e-3. 13,09©- 3, 11,44©- 3, 10,08©- 3,0 .94© -3, 

& 8. 0-3. 7. 20-3, 6. 53e-3, 5. 970-3,5.940-3/ 

data gams t/0 . 99987 , 0 . 999999999 . 0 . 99973 , 0 . 999 13,0. 99823 ,0 . 99708 , 

& 0.99568,0.99406.0.99225,0.99025,0,98807/ 

1 f ( to .gt . 50. ) go to 10 
Ita-lf ix(ta*.2)+1 
f rac“ta-f 1 oat (5*( 1 ta- 1 ) ) 

1 ta 1" i to+ 1 

8Ut"( sutt( 1 ta 1 ) -sutt( i ta) ) ♦0.2*f'rac+autt( 1 ta) 
nu*«(nut( ito1)-nut( i ta) ) *0, 2*f rac+nut( Ita) 
gamsY^«( (gamst( 1 tal )-gamst( 1 ta) ) ♦ . 2*f rac+gamst( 1 ta) )*980. 
return 

10 sut*»sutt( 1 1 ) 
nu*nut( 1 1 ) 
gamsw"gamst( 1 1 ) 
return 
end 
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DOCUMENTATION OF THE COMPUTER PROGRAM SPLASH 
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Documentation of the Computer Profiram SPLASH 

A complete documentation of the computer program SPLASH . FORTRAN is given 
by Mllly and Eagleson (1980). Here, only the procedure to achieve convergence 
of the results will be described and the way of attaching a file to it, in- 
cluding the boundary conditions of the area under investigation. 

Two parameters were varied in order to achieve convergence. Those were; 

1. XERR 

The parameter XERR represents the maximum allowed change of soil-moisture 
at every node and at every time-step. That is, 

■ nSes 1 

As this parameter decreases, the accuracy of calculations increases. In 
studying the catchments of Santa Paula and Clinton, it was found that 
convergence of the results occurs when XERR = 0.0005. 

2. ZRAT 

The parameter ZRAT is given by: 

_ (length of top element) x (number of elements) 

(total column length) 

For a fixed number of nodes, the value of ZR/vT was varied until satis- 
factory convergence was achieved. The results for Clinton, Massachusetts 
and Santa Paula, California are shown in Figures 31 and 32, respectively. 

It was found that satisfactory convergence is achieved for Clinton, 
when XERR = 0.0005, ZRAT = 0.01 and n = 21. For Santa Paula it was found 
that convergence can be considered achieved when XERR = 0.0005, 

ZRAT = 0.01, and n - 41. 
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Inputs for SPLASH 


First the program "input" described by Milly and Eagleson (1980) must be 
run, in which the number of nodes and the manner of setting up the nodes is 
established, and also the parameters XERR and initial are defined. 

By running "input" File 98 is created. This file must then be combined 
with a file including the soil properties and the atmospheric boundary condi- 
tions of the area under investigation. This file is also described with de- 
tails by Milly and Eagleson (1980). By combining those two files. File 15 
is created. 

Then, the program SPLASH . FORTRAN is ready to run, using as input File 
15. 
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